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ABSTRACT 

We perform local, vertically stratified shearing-box MHD simulations of protoplanetary disks (PPDs) 
at a fiducial radius of 1 AU that take into account the effects of both Ohmic resistivity and ambipolar 
diffusion (AD). The magnetic diffusion coefficients are evaluated self-consistently from a look-up table 
based on equilibrium chemistry. We first show that the inclusion of AD dramatically changes the 
conventional picture of layered accretion. Without net vertical magnetic field, the system evolves 
into a toroidal field dominated configuration with extremely weak turbulence in the far-UV ionization 
layer that is far too inefficient to drive rapid accretion. In the presence of a weak net vertical field 
(plasma (3 ^ 10^ at midplane), we find that the magnetorotational instability (MRI) is completely 
suppressed, resulting in a fully laminar flow throughout the vertical extent of the disk. A strong 
magnetocentrifugal wind is launched that efficiently carries away disk angular momentum and easily 
accounts for the observed accretion rate in PPDs. Moreover, under a physical disk wind geometry, all 
the accretion fiow proceeds through a strong current layer with thickness of ~ Q.iH that is offset from 
disk midplane with radial velocity of up to 0.4 times the sound speed. Both Ohmic resistivity and 
AD are essential for the suppression of the MRI and wind launching. The efficiency of wind transport 
increases with increasing net vertical magnetic flux and the penetration depth of the FUV ionization. 
Our laminar wind solution has important implications on planet formation and global evolution of 
PPDs. 

Subject headings: accretion, accretion disks — instabilities — magnetohydrodynamics — methods: 
numerical — planetary systems: protoplanetary disks — turbulence 



1. INTRODUCTION 

Protoplanetary disks (PPDs) are gaseous disks sur- 
rounding protostars. The gas in PPDs are found to be 
rapidly accreting to the protostar with accretion rate of 
\q-^±^Mq yr~^, with typical disk lifetime of about 1-10 
Myrs (e.g., iHartmann et al.l [l998t ISicilia-Aguilar et al] 
|2006[ ). Despite large number of observational programs 
aiming at revealin g the structure, composi tion and evolu- 
tion of PPDs (see IWilliams fc Ciezal[20Tll and references 
therein), two crucial theoretical questions on the gas dy- 
namics of PPDs remain poorly understood: What is the 
level of turbulence in PPDs? How efficient is angular 
momentum transport in PPDs? The answer to these 
questions arc essential to understanding the structure 
and evolution of the PPDs, as well as a series of pro- 
cesses in_2lajie^JonnaAi^ growth 
(e.g., iBirnstiel et all 120101: IZsom et al.l 120101) . tran sport 
of solids (e.g.. lGaraudll2007t iHughes fc Armitagell2010l ) 
are both sensitive to the radial structure of PPDs and 
level of turbulence. Current models for p lanetesimal for- 
mation such as the stream ing instability (jJohansen et al.l 
I2009| : |Bai fc Stond l2010allbl ) and gravitational instabil- 
ity (|Ybudin 2 0111) generallv favor weak turbulence and 
small radial pressure gradient. Moreover, turbulent mix- 
Electronic address: xbai@cfa.harvard.edu 

^ Department of Astrophysical Sciences, Peyton Hall, Princeton 
University, Princeton, NJ 08544 

Institute for Theory and Computation, Harvard-Smithsonian 
Center for Astrophysics, 60 Garden St.. MS-51, Cambridge, MA 
02138 

3 Hubble Fellow 



ing and dis k winds have a s ignificant influe nce on the dis k 
chemistry (jSemenov et al.l [2OO61 : Hcinzeller et al.l 1201 It) . 
When planets have formed, planet migration via planet- 
disk interaction also depends on the disk radial profile 
and diffusion proces ses (e.g., iPaardekooper et aITl2010l : 
iBaruteau et al.ll201lD . 

1.1. The Current Understanding of Accretion in PPDs 

We arc mainly interested in the T-Tauri (class II) phase 
of PPDs when the envelope infall has ended and the en- 
tire disk is visible. At this stage the disk is in general not 
massive enough f or gravitational instability to take place 
(|Zhu et al.ll2Cil0t ). In this paper, we focus on magnetic 
mechanisms. 

Most regions of PPDs'^ are ver y weakly ionized 
(|Havashil98ltlIgea fc Glassgoldlll999f ). hence the gas dy- 
namics is strongly affected by non-ideal magnetohydro- 
dynamics (MHD) effects due to the finite gas conductiv- 
ity, which include Ohmic resistivity. Hall effect and am- 
bipolar diffusion (A D). All three effects are r elevant and 
important in PPDs (|Wardlell2007HBaill2011a[) . Generally 
speaking, Ohmic resistivity dominates in dense regions 
with weak magnetic field (e.g., midplane in the inner re- 
gion of PPDs), AD dominates in tenuous regions with 
strong magnetic field (e.g., disk surface and outer region 
of PPDs), while the Hall regime lies in between. 

* The region that is close to the inner edge of PPDs is sufBciently 
hot (> lO'^K) due to direct illumination by the star that thermal 
ionization of Alkali species Na and K wil l provide sufficient ion- 
izatio n and the gas behave as ideal MHD ijUmebavashi fc Nakanol 
1 19881) , which is not the concern of this paper. 
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It is widely believed that PPDs arc turbulent as 
a result of the magn etorotational instability (MRI, 
iBalbus fc Hawlevlll998[) . The MRI turbulence transports 
angular momentum radially within the disk that allows 
the majority of the materials to be accreted onto the 
protostar while a small fraction of mass disperses away. 
Non-ideal MHD effects in PPDs strongly modify the be- 
havior of the MRI. Currently, most studies of the MRI 
in PPDs take into account only the effect of Ohmic re- 
sistivity, and it is found that in the inner region of PPDs 
(about 0.5-5 AU), the disk midplane is too we akly ion- 
ized f or the MRI to operate (i.e., the dead zone. iGammid 
Il996| ) , while the disk surface is still prone to the MRI and 
should be turbulent (i.e., the active layer). A large num- 
ber of numerical simulations have been conducted either 
in the local shearing-box framework or using full global 
approach to study and characterize the gas dynamics of 
the active layer and the dead zone, as well as exploring 
their physical cons equences (e.g. , Fleming & Stone 200^ 
Turner ct al."2007': 'Turner & Sano 2008; Ilgner & Nelson 
2008: Qishi fc Mac Low 2009: Dzyurkevich et al. 2010; 
Hirose fc Turneij 120111: IPlaig et al.l 120121: IPlock et al.l 



20121 ). These studies show that the boundary between 



the active layer and the dead zone is characterized by 
the Ohmic Elsasser number A « 1 (see equation ([6]) 
). Moreover, the dead zone is not completely "dead" 
in the sense that sound waves injected from the active 
layer bounce back and forth and give rise to some small 
Reynolds stress. Nevertheless, angular momentum trans- 
port is largely dominated by the MRI turbulence in the 
active layer, and the strength of the turbulence (consid- 
ering Ohmic resistivity only) appears to be able to drive 
rapid accretion consistent with observations. 

Comparatively, the effects of non-ideal MHD terms 
other than Ohmic resistivity on the MRI in PPDs are 
less well understood. They are extensively studie d in the 
linear regime without verti c al stra ti fication (see iWardlel 
19991: IBalbus fc Tergueml 120011: IWardle fc SalmerT^ 



2012 for the Hall effect. iBlaes fc BalbusI 119941: 



Kunz fc Balb^'200l 'Desch|[200jfor the effect of AD 



and [Pandey & Wardle 2012| for a general study ) , and 
with vertical stratification (jSalmeron fc Wardld l2005l 
[200l . In the non-linear regime, so far all numerical 
simulations (all using the local shearing-box approach) 
focus on individual non-ideal M HD e ffects such as MRI 
with Ohmic and Hall terms (ISano fc Stone 2002a 
and M RI with AD (jHawlev fc StoneHl998l: iBai fc Ston^ 
1201 It ). Most of the simulations are vertically unstrat- 
ified. These simulations provide useful criteria for the 
MRI to be self-sustained in the non-linear regime and 
the results were applied in the framework developed by 
iBail ()2011al) to estimate the efficiency of MRI-driven an- 
gular momentum transport in PPDs. It was shown that 
even in the most optimal scenario, the MRI-driven ac- 
cretion rate falls below typical observed rate by about 
an order of magnitude at the inner disk around 1 
AU. The main reason is that the strength of the MRI 
turbulence is expected to be dramatically reduced in 
the conventional "active layer" of the disk once AD 
is taken in to a ccount. Similar conclusions were also 
drawn from iPer cz-Becker & Chiang (2011 a,b); a nd from 
iMohantv et all (|2013D : IDzvurkevich et al.l (j2013t) for dif- 



ferent stellar masses. Relatively large accretion rate can 
be achieved in the outer region of the disk under optimal 
magnet ic field geo metry, and with the assistance of tiny 
grains (|Baill2011bD . 

An alternative scenario for describing the gas dy- 
namics in PPDs is the picture o f magnetocentrifu- 
gal w ind (jBlandford fc Pavni 119821 : iPudritz fc NormanI 
Il983t ): outflowing gas from accretion disks can be cen- 
trifugally accelerated along magnetic field lines when 
the inclination angle of the poloidal field is above 30° 
(relative to the disk normal)^. The magnetocentrifugal 
wind scenario has also been extensively explored with 
global simulations. Early global simulations treat the 
disk as a boundary condition (i.e., razor-thin) with ax- 
isym metry, and prescribe th e rate of outflow from the 
disk (jOuved fc Pudritzl 119971: iKrasnopolskv et al.l [l99l 
l2003t ). These simulations demonstrated the robustness 
of the magnetocentrifugal acceleration and collimation, 
and further found that the flow structure is sensitive to 
the prescribed rate of mass loading from the disk, and 
may lead to episodic form a tion of jets (iOuyed et al.|[l997l : 
IKrasnopolskv et al.l 119991: lAnderson et al.ll2005t) . More 
recent simulations (most of which are two-dimensional) 
that do resolve the disk generally rely on artificially pre- 
scribed and excessively large diffusion (that is unjusti- 
fied) to prevent rapid magnetic flux accumulation near 
the central objec t as mass a ccretes (|Kato et al.l 120021 : 
ICasse fc KeppensI l200t I2004D . and the resulting wind 
properties largely depend on the pr escribed resistivity 
profile in the disk (iZanni et al.ll2007D . 

In reality, the wind launching process is governed by 
the microphysics within the disk, with mass loading 
rate determined by requiring that the flow smoothly 
passe s t he slow magnetosonic point (I Wardle fc Koenigll 
Il99l 119951 : lOgilvie fc Li^lMlt lo'gilviell2012D . It 
was found that launching a laminar disk wind gener- 
ally requires a relatively strong vertical background mag- 
netic fleld with around equipartit i on strength at disk 
midp lane (jWardle fc Koenigll 119931 : iFerreira fc Pelletierl 
|1995[ ). though weaker fleld is possible in the presence of 
strong ambipolar diffusion (Li 1996). The vertical mag- 
netic field can not be too strong, which would make it 
difficult to be bent by the disk material, and would re- 
sult in substanti al sub-Keplerian rotation that hinders 
wind-launching (jShu et al.l 120081: lOgilvid l2012l) . More 
detailed study of the wind- launchin g criteria and repre - 
sent ative solutions were pr esented in lKonigl et al.l ()2010l ) 
and lSalmeron et aD (j2011t ) where all the non-ideal MHD 
effects were taken into account. It appears that the 
MRI and the lamin ar wind scenarios tend to mutually 
exclude each other (jSalmeron et al.l I2007D : MRI oper- 
ates when the vertical magnetic field strength is well be- 
low equipartition strength at the disk midplane, while 
launching an magnetocentrifugal wind requires the verti- 
cal field strength to be around equipartition. The strong 
magnetic field in the magnetocentrifugal wind scenario 
would drive very efficient accretion, making the disk 



5 The X-wind model llShu et aLlligQ^ I is also magnetocentrifugal 
in nature, with the wind launched near the inner edge of the disk. 
We are interested in the wind launched from radially extended 
region in PPDs. 
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much more tenuous (with very small surface density) 
than that in the MRI scenario for a standard accretion 
disk (jCombet fc Ferreirall2008[) . 

To briefly summarize, MRI generally requires the back- 
ground vertical field to be weak, and has difficulty in ac- 
counting for the rapid accretion rate in PPDs, especially 
in the inner region of ^ lAU. The alternative picture 
of magnetocentrifugal wind generally requires the pres- 
ence of strong background vertical field of equipartition 
strength that either drives accretion that is too rapid or 
results in a tenuous disk, but the microphysics of the 
wind launching process still requires more realistic treat- 
ment. 

1.2. This Work 

We conduct vertically stratified shearing-box simula- 
tions of a local patch of a PPD that include a realis- 
tic profile of both Ohmic resistivity and AD coefficients. 
The diffusion coefficients are interpolated from a pre- 
computcd look-up table in real simulation time based 
on the gas density, temperature (fixed) and ionization 
rate (calculated from the density profile). We have not 
included the Hall effect (which is numerically more chal- 
lenging and demanding), and this is the first step to- 
wards understanding the non-ideal MHD of PPDs be- 
yond Ohmic resistivity. This is the first time that Ohmic 
resistivity and AD are simultaneously included in ver- 
tically stratified simulations to study the gas dynamics 
of PPDs that incorporates the disk microphysics in the 
most realistic manner. 

All our simulations except one include a weak verti- 
cal magnetic field, which is likely to be more realistic for 
a local patch of a PPD. The other reason for including 
a vertical field is that such a field geometry is more fa- 
vorable for the MR I to operate in the presence of AD 
()Bai &: Stondl20l"l[ ). We note that most vertically strat- 
ified shearing-box simulations to date adop t a zero net 
vertical magnetic flux field geoinetry (e.g., IStone et"al] 
[1991 iShi et aLllMot iDavis et al.|[2oTol ) and we will show 
that such field geometry would make the strength of the 
MRI turbulence diminishingly small due to AD. Includ- 
ing net vertical magnetic flux in shearing-box simulations 
places strong demands on the robustness of numerical 
algorith ms, especially in th e magnetic dominated disk 
corona (jMiller fc Stonel[2000( ) . Recently such simulations 
are successfully performed by several groups, mostly in 
the ideal MHD regime, and it was shown that the in- 
clusion of net vertical magnetic field always leads to an 
outflow from an MRI-turbulent disk ( Suzuki fc Inutsukal 
[2009t iFromang et all [20T1 iBai fc Stone) 12011) . In the 
context of PPDs, it was found that including a net ver- 
tical magnetic flux does not change the basic picture of 
layered accretion (again, considering Ohmic resistivity 
only), but stronger net vertical flux leads to stronger 
MRI turbulence in the active layer and reduces the ver- 
tical extent of the dead zo ne, as well as stronger outflow 
(jOkuzumi fc Hiros^l2011f ). Although these results sug- 
gest the simultane ous existence of the MRI and magne- 
tocentrifugal wind, iBai fc Stond (|2013[ ) pointed that the 
outflow from MRI active region is unlikely to be directly 
connected to a global magnetocentrifugal wind due to the 
MRI dynamo and symmetry considerations. Our simu- 



lations in this paper will further address the potential 
connection of the disk outflow to an magnetocentrifugal 
wind in the context of PPDs, and we arrive at positive 
co nclusions, in contras t with the ideal MHD case studied 
bv lBai fc Sto "^ ([20l3h . 

More specifically, we consider a fiducial model that 
corresponds to a minimum-mass solar nebular at 1 AU, 
assuming solar abundance of chemical composition and 
0.01% (in mass) of 0.1/im sized grains for the chemistry 
calculation. We find that although the initial condition 
is unstable to the MRI (with net vertical magnetic field 
much weaker than equipartition), the disk rapidly ad- 
justs to a new laminar configuration that is stable to 
the MRI. The new laminar state is characterized by an 
outflow launched by the magnetocentrifugal mechanism, 
and the outflow can achieve a physical wind geometry 
(poloidal streamlines at the top and bottom sides of the 
disk bend towards the same radial direction) by having 
a thin current layer where the horizontal field flips. The 
magnetocentrifugal wind launched in this scenario can ef- 
ficiently transport angular momentum to account for the 
observed PPD accretion rate while without being too ef- 
ficient to deplete the disk (as in the conventional wind 
scenario). For clarity, we focus on the physical properties 
of the new laminar solution and together with a param- 
eter study all at 1 AU in this paper, extension of the 
results to other disk radii will be presented in a compan- 
ion paper. 

This paper is organized as follows. In Section [21 we 
describe the numerical method and the overall setup of 
our simulations. Three contrasting simulations are pre- 
sented in Section |3l highlighting the new laminar wind 
solution. In Section 21 we study the physical properties 
of the new laminar wind solution in detail and discuss the 
wind launching mechanism. In Section [SI we conduct a 
thorough parameter study to further explore the proper- 
ties of the laminar wind. We discuss the robustness and 
implications of our wind solution in Section [51 together 
with the conclusion. 

2. SIMULATION SETUP 

2.1. Formulation 

We use the Athena MHD code ([Stone et al.ll2008D to 
study the gas dynamics in PPDs and perform three- 
dimensional numerical simulations. We consider a local 
patch of a PPD and adopt the conventiona l shearing-box 
approach ([Goldreich fc Lvnden-Beilll965[ ). MHD equa- 
tions are written in a Cartesian coordinate system in the 
corotating frame at a fiducial radius with Keplerian fre- 
quency Q,: 

(pu) = , (1) 



dpu 



dt 



' V • {pu^u + T) 



2m X n -I- xe^ - VL ze^ 



inhere T is the total stress tensor 
T= {P + B^/2) I- 



B^B 



(2) 
(3) 



I is the identity tensor, p, it, P are the gas density, 
velocity and pressure respectively, B is the magnetic 
field, ex,ey,ez are unit vectors pointing to the radial. 
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azimuthal and vertical directions respectively, where $7 
is along the direction. Note that the equations are 
written in units such that magnetic permeability is 1. 
Vertical gravity is included to account for density strati- 
fication. Periodic boundary conditions are used in the 
azimuthal direction, while the radial boundary condi- 
tions are shearing periodic as usual. Shearing-box source 
terms (Coriolis force and tidal gravities) have bee i i read - 
ily implemented in Athena ( Stone fc Gardiner! 120100 . 
which uses an orbital advection scheme that splits the 
system into an advective part for the background shear 
flow — 3f2x/2ej,, and a fluctuating part with velocity fluc- 
tuation v. 

3 

V = M + -fixCj, . (4) 

The advection scheme not only accelerates the calcula- 
tion by permitting larger time steps, but also improves 
the accuracy by making the truncation error independent 
of radial location. 

We use an isothermal equation of state P ~ pc^, where 
Cg is the isothermal sound speed. In reality, multiple ra- 
diative processes such irradiation and scattering are in- 
volved in PPDs, m aking the disk surface layer substan- 
tially hotter (e.g., iHirose fc funi er 20 111). The hotter 
disk surface affects the properties disk wind, however, 
the main problems addressed in this paper are largely 
magnetic: the suppression/survival of the MRI and wind 
launching arc largely controlled by non-ideal MHD ef- 
fects, where thermodynamics only plays a minor role. 

Being weakly ionized, PPDs are not perfectly conduct- 
ing, which is reflected in the non-ideal MHD terms in the 
induction equation (e.g.. IWardlil2007t lBaill2011al ) 

dB 

— = Vx{uxB)-Vx[r]oJ+r]H{J^B)+7uJ±] , (5) 

where J = V x S is the current density, B denotes 
unit vector along B, subscript "j^" denotes the vector 
component that is perpendicular to B, r]o,ilH and rjA 
are the Ohmic, Hall and the ambipolar diffusivities. Note 
that u represents the velocity for the bulk of the gas 
(i.e., neutrals), while tracer amount of charged species 
provides conductivity and gives rise to non-ideal MHD 
effects. 

The magnetic diffusivities depend on the number den- 
sity of the charged species, and are characterized by the 
dimensionless Elsasser numbers, defined as 



A = 



Ha = 



— "A 



At 



— ""A 



VA^ 



(6) 



for Ohmic, Hall and AD respectively, where va — \/B^/p 
is the Alfven velocity. In the absence of grains, we have 
A oc B^, Ha oc B and Am being independent of B (see 
next subsection for details). Generally speaking, self- 
sustained MRI turbulence requires these Elsasser num- 
bers to be greater than 1. 

Ohmic resistivity and ambipola r diffusion (AD) 
have been implem ented in Athena ([Davis et al.l 120101 : 
iBai fc Ston^ 1201 It ). Furthermore, we have imple- 
mented super time-stepping that substantially acceler- 
at es the calcul a tion, which we have described in detail 
in lSimon et al.l (|2013[ ). Although our simulations do not 



include the Hall term, we do evaluate r]H and assess its 

importance in our analysis. 
We use natural unit in our simulations, where Cs = 1, 
= 1. The initial density profile is taken to be Gaussian 



poexp(-;2V2H2) 



(7) 



with po = 1 being the midplane gas density in code unit, 
and H = Cs/il = 1 is the thermal scale height. We 
perform simulations with both zero and non-zero net- 
vertical magnetic flux. For simulations with zero net 
vertical flux, the initial field is purely vertical given by 
Bz = i?! sin 27ra;/i/a; , where is the radial size of the 
simulation box, and Bi parameterized by the midplane 
plasma /3i = 2poc^/Bf. For simulations with net-vertical 
flux, we add a uniform vertical field Bq (parameterized 
by midplane plasma /3o = '^Poc^/ Bq) on top of the sinu- 
soidally varying component. 

Physically, we consider a patch of the PPD using the 
minimum-mass solar nebular (MMSN, IWeiden schilliii"^ 
Il977tlllavashilll98l[ l disk model, with the surface density 

given by S(i?) = 1700i?^u^^g cm~^, and temperature 

— 1/2 

given by r = 280i?^u K, where Rav is the radius to 
the central star (whose mass is fixed at IMq) measured 
in astronomical unit (AU). The mean molecular weight 
of the neutrals is taken to be pn = 2.34m/f from which 
the sound speed Cg ~ \/ kT / pnmH = 1 km s^^ (at 1 
AU) hence other quantities can be easily evaluated. In 
particular, i3 = 1 in our code unit corresponds to field 
strength of 13.1 Gauss. These physical scales are needed 
to normalize the magnetic diffusivities (next subsection) 
to code unit. 

We use outflow boundary condition in the vertical di- 
rection which copies the density, velocity and magnetic 
fields in the boundary cells to the ghost zones, with the 
density attenuated following the Gaussian profile to ac- 
count for vertical gravity. In the case of mass inflow, 
the vertical velocity is set to zero at the ghost zones. A 
density floor of ppioor = 10~^ (in code unit) is applied to 
avoid numerical difficulties at magnetic dominated (low 
plasma /3) regions. We have checked that horizontally 
averaged densities in the saturated states of all our sim- 
ulations are always well above the density floor^. More- 
over, the use of outflow boundary condition no longer 
conserves mass in the simulation box. We compensate 
for the mass loss so that a stea dy state can be achi eved, 
following the same procedure as lBai fc Stond (|2013( ): the 
density of each cell is modified by the same proportion 
at the end of each time step so that the total mass in 
the simulation box remains the same. For most of our 
runs, the mass change over the duration of our simula- 
tions (if mass conservation were not enforced) is only a 
tiny fraction of the total mass. 

2.2. Calculation of Magnetic Diffusivities 

The magnetic diffusivities are evaluated based on the 
chemistry calculation. Instead of evolving the chemi- 
cal network i n real time, as don e by a number of pre- 
vious works (jTurner et al.l l2007t iTurner fc Sanol 120081 : 

^ Except for Run OA-nb-F3 where a density fl oor o f ppioor = 
10~* is applied, which will be discussed in Section 13.21 
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llkner fc NelsM][2008h. we assume equilibrium chemistry 
(similar to iHirose &: Turneil l2011| ). because the recom- 
bination time has been s hown to be much shorter than 
the dynamical time scale ()Baill2011al). We adopt a com- 
plex chemical reaction network signer fc NelsonI [20061 : 
iBai fc Good man! 120091: lBaill2011arthat is based on the 
UMIST database IWoodall et al.ll2007t) . In our fiducial 
model considered in this paper, we fix the elemental com- 
position to be solar, with well-mixed 0.1/im grains whose 
abundance is = 10~^ in mass (about 0.01 solar, corre- 
sponding to substantial grain growth and settling). We 
follow the sa nie procedur e and methodology described 
Section 3.4 in iBail ()2011a[ ) to evolve the network for 10^ 
ye ars, with further detail s provided in Sections 3.2-3.5 
of iBai fc GoodmanI ()20Q9l ). Given the chemical compo- 
sition, variable parameters of the network include gas 
density p, gas temperature T, and the ionization rate 
^, which are scanned to give a complete coverage of the 
parameter space relevant to PPDs. The outcome of the 
scan is a look-up table of magnetic difFusivitics that is 
read into the code so that r/o, fin and rjA at each grid 
cell can be evaluated by interpolation in real simulation 
time. Since we adopt an isothermal equation of state, T 
is fixed, our look-up table is essentially two-dimensional 
{P and C).^ 

The ionization rate in the disk depends on the column 
density to the disk surface. We calculate the horizon- 
tally averaged vertical density profile in real simulation 
time, from which a column density profile can be re- 
constructed, an approach similar to previous works (e.g. 
[Turner et al . 2007). The sources of ionization include ra- 
dioactive decay, cosmic ray and stellar X-ra y ion i zations , 
with prescriptions given in Section 3.2 of iBail ()2011al ). 
with fixed X-ray luminosity of 10'^° erg s~^ and X-ray 
temperature of 5keV. In addition, we consider the ef- 
fect of far-ultraviolet ( F UV) io nization. According to 
iPerez-Becker fc ChiangI (j2011bf ). FUV photons almost 
completely ionize tracer species such as C and S and give 
ionization fraction of the order / = 10^^ — 10^"' with pen- 
etration depth 0.01 — 0.1 g cm~^ depending on the effec- 
tiveness of dust attenuation and self-shielding ^. For sim- 
plicity, we assume an ionization fraction of / = 2 x 10"^ 
in the form of carbon in the FUV layer, whose column 
density is chosen by default as Spuv = 0.03 g cm^^. 
Within the FUV layer, the magnetic diffusivities ex- 
pressed in the form of Elsasser numbers under the MMSN 

There are large uncertainties associated with the FUV ioniza- 
tion. For example, it has been noted that the FUV photons may 
be shielded if a dusty win d is launched from the inner disk near the 
dust sublimation radius l IBans &: Koniglll2012l ). Moreover, photo- 
chemistry in the FUV layer plays an important role in determining 
the m olecular composition and and ion abundances llWalsh et al.l 
120121) . and particularly, the FUV penetratio n depth may be over- 
estimated in lPerez-Becker fc ChiaM l|2011bl ') due to the conversion 
of C"*" and into molecular ions which facilitates recombination 
IIAdamkovics et al.l 120111 . Al Glassgold, private communication). 
Our treatment should be considered as a first approximation that 
captures basic physics rather than deep into the details. 



disk model arc found to be 



Am 



3.6 X 10' 
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- 1 -^Au^"^ 



Pa J 



Ha = 



6.4 X 10' 



10- 



-.R 



-1/8 
AU 



(8) 



where /3mid is ratio of magnetic pressure to midplane gas 
pressure. A smooth transition (across about 4 grid cells) 
on the magnetic diffusivities from the FUV ionization 
layer at S < Spuv to the X-ray/cosmic-ray dominated 
ionization layer (based on our chemistry calculations) at 
E > Spuv is applied. As Ohmic resistivity plays essen- 
tially no role in the low density region of the FUV layer, 
we simply do not reset Ohmic resistivity in this layer. 

We calculate the magnetic diffusivities from the num- 
ber density of all charged species at the end of the chem - 
ical evolution following Section 2 and 3.5 of lBail (j2011al ). 
Note that the Ohmic resistivity rjo is independent of 
magnetic field strength B, while Hall and ambipolar dif- 
fusivities do depend on B. In the grain-free case, we 
have TjH oc B, and tja oc B^. In this case, we can 
simply fit the proportional coefficients, Qh and Qa re- 
spectively, and put them into the look-up table. In the 
presence of small grains, a situation studied in detail in 
iBail ()2011b[ ). rjH (va) is proportional to B (B^) when 
B is sufficiently weak or sufficiently strong respectively, 
while is roughly proportional to B^^ {B'^) at some in- 
termediate field strength. In this case, we include in the 
look-up table the two proportional coefficients Qhii Qh2 
{Qai, Qa2) at weak and strong field regimes from the fit- 
ting respectively, together with a transition field strength 
Bi so that 



QhiB , 

QH2BfB~^ 



;h2 



B 



B < y/ QH2/QHlBi , 

VQh2/QhiB, <B<B, , (9) 
B>B, , 



VA 



'QaiB^ 

Qa2b: 

.Qa2B^ 



2 

i ' 



B < ^Qa2/ QaiB^ , 
^Qa2/QaiB^ <B<B, , (10) 
B>B, , 



By comparing with Figure 1 oflil (I2011b[ ). we see that 
the transition field strength Bi corresponds to the sit- 
uation that the ion gyro-frequency equals its collision 
frequency with the neutrals (or the ion Hall parameter 
equals one), which can be directly calculated given the 
gas density. 

Being diffusive processes, large Ohmic resistivity near 
the disk midplane and strong AD in the tenuous disk 
corona would significantly limit the code efficiency even 
super time-stepping is used to accelerate the calculation. 
In the Ohmic regime, the diffusive time step scales as 
b? lr\o where A is the minimum grid spacing. In prac- 
tice, we add a cap to the Ohmic resistivity so that in code 
unit 770 < 10. Since the magnetic field strength near the 
disk midplane never reaches equipartition in all our sim- 
ulations, the Elsasser number at the disk midplane is 
always smaller than 0.1, well below the threshold value 
1. Also, this cap value of 770 makes the diffusion time 
scale much smaller than the dynamical time scale, hence 
captures the basic effect of strong diffusion at disk mid- 
plane even if resistivity is much higher in reality. Note 
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that our resistivity cap is much larger than most pre- 
vious works (thanks to the use of super time-stepping), 
where the cap was of the order 0.01 in natural un it (e.g., 
[Fleming fc Stonell2003l: [Okuzumi fc Hirosell201ir ). In the 
AD regime, the time step scales as A'^-Am-/3 ( /3 is the lo- 
cal ratio of gas to magnetic pressure). In the same spirit, 
we apply a floor to Am so that Am • /3 > 0.1 in every grid 
cell. This floor value of Am is again sufficiently small so 
that it does not make the otherwise stable field configu- 
ratio n unstable to the MRI (sec Figure 16 of lBai fc Strn^ 
1201 It) , and it retains the effect of strong diffusion. 

3. FIDUCIAL SIMULATIONS AND RESULTS 

In this section, we present three benchmark simula- 
tions with very similar initial setup but evolve into dra- 
matically different states. All three simulations are three- 
dimensional (3D) shearing-box with vertical stratifica- 
tions, located at 1 AU in a MMSN disk, adopting a chem- 
istry model with well- mixed OAfim grains with abun- 
dance of 10~^ (1% solar), and a FUV column density 
^Fuv = 0.03 g cm~^. All simulations arc run for about 
150 orbits (900rj-i). 

In the first simulation (Run 0-b5), only Ohmic resistiv- 
ity is included, which aims at modeling the conventional 
picture of layered accretion. We have also included a 
net vertical magnetic flux of (3o = 10^. A sinusoidally 
varying vertical field of /3i = /3o/4 is also added as initial 
condition to avoid the strong channel flows (but has no 
effect on the saturated state of the system). The sim- 
ulation box size is AH x 8H x 16H in the radial (x), 
azimuthal (y) and vertical (z) dimensions respectively, 
with a computational grid of 96 x 96 x 384 cells. Our 
simulations have relatively high resolution in x and z (24 
cells per H, or 34 cells if one defines the scale height to 
be s/2H, as in a number of works) to properly resolve the 
MRI turbulence (iDavis et al.ll2010l : ISorathia et al.ll2012l : 
iBai fc Stonell2011f ). and have relatively large simulation 
box to c apture the mesosca le structures of the MRI tur- 
bulence (jSimon et al.ll2012D . 

In the second simulation (Run OA-nobz), both Ohmic 
resistivity and AD are included, and we adopt a zero 
net vertical magnetic flux configuration with /3i = 1600. 
Again, the choice of f3i has no effect on the saturated 
state of the system. After experimenting with several 
initial setups, we find that due to the extremely weak 
level of turbulence, the gas near the disk surface has very 
limited magnetic support hence drops very rapidly, and 
we must reduce the box height to 12H (and 288 cells) 
so that the gas density at vertical boundaries is not too 
small which would severely reduce the numerical time 
step. Furthermore, we started the simulation with den- 
sity fioor pfloor — 10~^, but then found that the density 
at vertical boundaries is always at the level of the density 
floor. We thereby keep reducing the density fioor until 
Pfloor = 10~* when we no longer find artificial features 
near the vertical boundaries. 

Finally, we conduct a contrasting simulation from the 
above two cases (Run 0A-b5). The initial setup of the 
simulation is exactly the same as in the first simulation 
(Run 0-b5) except that AD is included. 

Figured] illustrates the initial profile of the Ohmic, Hall 
and ambipolar Elsasser numbers for our runs 0-b5 and 



TABLE 1 

Summary of All Fiducial Simulations. 



Run 


Diffusion 


/3o 


Box Size 


Section 


0-b5 


Ohm 


10^ 


4H x8H X 16H 


3.1 


OA-nobz 


Ohm, AD 


oo 


iH x8H X 12H 


3.2 


OA-b5 


Ohm,AD 


10^ 


iH x8H X 16H 


3.3, 4 


S-OA-b5 


Ohm,AD 


10^ 


H/8 X H/A X WH 


4.4 


E-OA-b5 


Ohm,AD 


10^ 


iH x8H X 8H 


4.4.1 


S-OA-b5-12H 


Ohm.AD 


10^ 


H/8 X H/4 X 12H 


4.5 


S-OA-b5-20H 


Ohm.AD 


10^ 


H/8 X H/4 X 20H 


4.5 


S-OA-b5-24H 


Ohm,AD 




H/8 X H/A X 2iH 


4.5 


MMSN disk model at lAU, X-ray luminosity of Lx = 10'"' ergs 
s~^ and temperature Tx = 5keV, well mixed 0.1/xm grains with 
abundance of egr = 10~*, penetration depth of 0.03g cm~^ for the 



FUV, ionization are, assumed for all runs . 

OA- bo (similar but not exactly the same tor run OA- 
nobz due to the different magnetic configuration). Also 
shown is the initial profile of plasma /3. From the disk 
midplane to surface, the dominant non-ideal MHD effects 
arc Ohmic resistivity. Hall effect and AD respectively for 
the initial field configuration, and MRI unstable region is 
located at around z = ±.2>H where all Elsasser numbers 
arc greater than 1 and plasma /3 is well above 1. The 
initial FUV ionization front is located at about z = zLAH. 
The density floor of pfioor = 10^^ is applied to regions 
beyond z = ±5i?, hence the Am and /3 curves flattens 
out (this artifact will disappear as the system evolves). 

The most common accretion diagnostics is the i?0 com- 
ponent of the Reynolds and Maxwell stresses, which mea- 
sures the local rate of radial transport of angular momen- 
tum 

Tr^ = + Tr^ = 7^ - B.By , (11) 

where the over bar indicates horizontal averaging. The 
total rate of radial angular m omentum transport is char- 
acterized by the a parameter (jShakura fc Sunvaevlll973l ) 

pdz 

In Figure [2] we show the time evolution of the hori- 
zontally averaged Maxwell stress —BxBy for the three 
fiducial simulations. Albeit for similar initial setup, the 
three runs show very distinctive characteristics as the 
systems evolve into saturated/steady states. The results 
from the three simulations will be analyzed in detail in 
the following three subsections. 

For clarity, we further provide a list of all our fiducial 
simulation runs and their parameters in Table [T] In par- 
ticular, we introduce the letter "S" for runs with very 
small horizontal domain, where the simulations are es- 
sential one-dimensional, and letter "E" for runs with en- 
forced even-0 symmetry (see Section^]). We use the label 
"bn" to denote plasma /3o — 10" for the vertical back- 
ground field. All other simulations are run for about 200 
orbits (1200O-1). 

3.1. The Ohmic- Resistivity- Only Run 

Run 0-b5 quickly develops into turbulence. From the 
upper panel of Figure [5J the separation between the 
highly turbulent active layer and the more or less quies- 
cent midplane region is clearly seen. We further show the 
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Fig. 1. — Initial Elsasscr number profile for Ohmic resistivity 
(red solid), Hall diffusivity (green dashed) and ambipolar diffusiv- 
ity (blue dash-dotted) for our fiducial runs 0-b5 and OA-b5 . Note 
the eap on Am and A around the disk midplane and on Am beyond 
about ±5H. Also shown is the initial profile for plasma /3 (black 
solid). Note that we show the Hall Elsasser number (Ha) while it 
is not included in the simulations. 

time and horizontally averaged vertical profiles of den- 
sity, magnetic pressure, Maxwell and Reynolds stresses 
in Figure [S] The time averages are performed from 
fit = 450 onward. 

The MRI turbulence generates strong magnetic field 
that buoyantly rises, and the disk surface quickly forms 
a strongly magnetical ly dominated corona beyond ±2H 
([Miller fc Stond [2000( ). as shown in the bottom panel 
of Figure |31 The gas density in the corona follows an 
exponential profile due to strong magnetic support, in- 
stead decreasing with height as an Gaussian in the hy- 
drostatic case. The velocity in the coronal regions is 
highly supersonic, with turbulent kinetic energy exceed- 
ing the gas pressure beyond ±5H. A strong outflow is 
launched f rom the active laye r of the disk as has been 
studied by ISuzuki et al.l (120101 ) , who also found that the 
mass outflow rate scales roughly linearly with the vertical 
net magnetic flux. In our simulation, the time and hor- 
izontally averaged mass outflow rate pv^ is found to be 
about 3.7 X 10~'^poCs from each side of the box, corn para- 
ble with the measurements in ISuzuki et al.l ()2010l) with 
the same (Sq. We note that the mass outflo w rate is not a 
well-c haracterized quantity in shearing-box ()Bai &: Stond 
l2013f ). these measurements should only be taken as a ref- 
erence. 

The disk midplane is too resistive to become MRI ac- 
tive, with the boundary of the active layer well described 



by the Elsasser number criterion: 



where vaz is the vertical component of the Alfven ve- 
locity. The Maxwell stress remains very small at the 
midplane for the flrst 20 orbits. However, the mid- 
plane magnetic field is then gradually amplified, while 
the flow in the midplane remains more or less lami- 
nar. We note that the Maxwell stress in the midplane is 
even larger than most regions in the active layer, which 
contrasts with some previous simulations with a similar 
setup, where t he Maxwell stress becomes very small near 
the midplane (jSuzuki et al.ll2010l: iHirose fc Turneij|2011l: 



lOkuzumi fc Hirosd [20Tll ) . The main reason for the dif- 
ference, while somewhat counterintuitive, lies in the us- 
age of a much larger resistivity cap in our simulations^. 
The large resistivity at the disk midplane in our simu- 
lations strongly suppresses electric current, leaving the 
horizontal magnetic field to be almost constant across 
the midplane (see the bottom panel of Figure ^ . There- 
fore, the midplane field strength is largely set by the field 
strength in the active layer. Note that although this is 
phenome nologically similar to the "undead" zone pro- 
posed by [Turner fc Sand ()2008() , it is conceptually very 
different. We note that it is essential for the resistivity 
cap to be large enough so that the gas and magnetic field 
become decoupled at the cap value, and we have further 
tested that the simulation results are independent of the 
resistivity eap once its value is greater than 3 or so. 

By contrast, the Reynolds stress pVxVy as well as ki- 
netic energy clearly have a large dip in the undead re- 
gion between z = ±2iJ. There is still random mo- 
tion near the midplane due to th e sound waves launched 
from the base of the act ive layer (jFleming fc Ston^l2003l : 
lOishi fc Mac Lowll2009() . though the velocity amplitude 
is at least one order of magnitude smaller than that in 
the active layer. 

Integrating the profiles of the Maxwell and Reynolds 
stresses using ([T^. we obtain the Shakura-Sunyaev pa- 
rameter aMax ~ 1-3 X 10~^, and aRcy ~ 1-2 x 10""^. 
Assuming steady state accretion, one can express the ac- 
cretion rate for a MMSN disk model as 



27raEc? 



i.2 X IQ-^aRj^l/^AlQ yr^^ 



(13) 



Therefore, in the absence of AD, and at the fiducial 
location i? = 1 AU, we obtain an accretion rate of 
1.2 X 10~^Mq yr~^, large enough to account for the 
observed accretion rates in most T-Tauri stars. 

Given the usage of more realistic resistivity profiles and 
chemistry models, as well as the much larger resistivity 
cap enabled by the super time-stepping technique, our 
run with pure Ohmic-resistivity deserves more discussion 
on its own right. Nevertheless, we only consider our run 
0-b5 as a test and reference case since our main point of 
interest is the effect of AD on the gas dynamics of PPDs, 
which makes a dramatic difference from the conventional 
picture of PPD accretion. 

3.2. Zero Net Vertical Flux Run with both Ohmic 
Resistivity and AD 

Setting the initial density floor of pfloor = 10~^, we find 
that the MRI turbulence sets in at the beginning, then 
the simulation gets stuck with relatively strong magnetic 
field (dominantly toroidal) accumulating near the verti- 
cal boundaries with relatively little activities. The fiux 
does not escape and provides very little magnetic support 
due to the flat proflle, and the density near the vertical 
boundaries stays at the value of pfloor- This situation 
implies that the gas near vertical boundaries tends to 
fall back towards disk midplane, but this is prevented by 

* See last paragraph of Section 2. In fact, we start the simulation 
with a smaller resistivity eap of rjo = 1 in code unit. At t = 
120Q~^, the eap is raised to its final value t]q = 10, and the 
Maxwell stress at disk midplane rises shortly afterwards. 
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Fig. 2. — The space-time plot for the (horizontally averaged) vertical profiles of the Maxwell stress —B^By in the three fiducial runs 
0-b5 (top), OA-nobz (middle) and OA-b5 (bottom). The colors are in logarithmic scales (in code unit). White contours in the upper panel 
correspond to vertical Elsasser number = Vj^^/rjo^ = 1- The discontinuous transitions in the middle panel correspond to the sudden 
reduction of the density floor which attains the final value of Pfioor = 10~* at i = 230f2~^. 




Fig. 3. — Vertical profiles of the Maxwell stress and Reynolds 
stress (upper panel), as well as gas/magnetic pressure and kinetic 
energy (lower panel) in the fiducial run with only Ohmic resistivity 
(0-b5). 

the outflow boundary condition. Therefore, we gradu- 
ally lower the density floor over the course of the simula- 
tion and find that the phenomena of artificial magnetic 
flux accumulation and density cutoff at pfioor disappears 
when pfloor is brought down to 10~^, which is achieved 



Fig. 4. — Same as Figure [3] but for zero net vertical flux run 
with both Ohmic resistivity and AD (OA-nobz). Note the different 
scales. 

at t = 230r2~"'^. The whole process is reflected in Figure 
[21 From this time onward, the system evolves smoothly 
and no artifacts arc seen from the vertical boundaries. 
Also note that setting such small density floor of 10^® at 
the beginning would lead to dramatically small simula- 
tion timcstep which is not quite realistic, hence gradual 
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reduction of pfloor is necessary, although one has to be 
patient to get rid of the artifacts. 

We see from Figure [5] that in the saturated state, the 
disk only has extremely weak level of MRI turbulence in 
a very thin layer between about z ~ 4 — 5H. Note that 
the FUV ionization front is located &i z — ±4iJ, beyond 
which we find Am > 50 and the gas behaves almost as in 
the ideal MHD regime (since density profile is still largely 
Gaussian, one can still use Figure [T] as a reference for the 
profile of Am), hence the MRI operates in this region. 
Below z — AH, the value of Am is of order unity or 
smaller, and there is no evidence of turbulent activity. 
This is consi s tent w ith unstratified simulation results of 
iBai fc Stoiii (|2011[ ). where it was found that the MRI 
is suppressed for Am < 1 in the absence of net vertical 
magnetic flux. 

There is some residual Maxwell stress around the disk 
midplane as a result of initial conditions that slowly di- 
minishes over the course of the simulation. We extract 
the profiles of various physical quantities and show them 
in Figure m where the time average is performed between 
fit =- 750 and 900. The MRI-active region exhibits as 
bumps in the stress plot, as well as the bumps in the 
kinetic energy plot. The magnetic field strength roughly 
stays constant through the entire disk, which is likely set 
by the turbulent activities in the active zone: the value 
of plasma /3 (ratio of gas pressure to magnetic pressure) 
reaches order unity within the turbulent layer. Moreover, 
we find that the magnetic field is predominantly toroidal 
in the entire disk. 

From the simulation, we find the Shakura-Sunyaev 
parameter using (fT2)) to be ofMax ~ 1-7 x 10~^, and 
aRey ~ 1-3 X 10^^. For such small level of stress, 
we find that the resulting accretion rate is only about 

2.5 X 10"^"'^ Mq yr"-'^, which is three orders of magnitude 
too small compared with observations. 

We can compare this result with optimistic predic- 
tions of the MRI-driven acc r etion r ate using the semi- 
analytical framework of iBail (|2011al ). Using this frame- 
work, we first extract the density profiles and the pro- 
files of Am and r\o- Assuming constant magnetic field 
strength across the MRI-active layer, the MRI-drive n ac- 
cretion rate can be expressed as (Equation (28) of iBail 
I2011ar ) 

81^2 Active a ' ^ ^ 

where the integral is performed in regions where t he MR I 
is permitted based on the criteria (20) of iBail ()2011al ). 
We scan the field strength to maximize M, which gives 

6.6 X WT^'^Mq yr~^ as an upper limit of the MRI-driven 
accretion rate. We see that this value is a factor of more 
than 20 larger than our simulation result. The main 
reason that our simulation yields an accretion rate that is 
significantly smaller than theoretical expectations is that 
the fi eld geometry is not optimal: in both the ideal MHD 
(e.g., iHawlev et all [l995l: IBai fc Stond [20ll and non - 
ideal MHD ('e.g.. lFleming et al.ll2000HBai h Stonell201lD . 
one finds stronger turbulence in the presence of a net 
vertical magnetic flux, and the optimistic estimate can 
possibly be achieved only when a optimal field geometry 
is realized. 



This simulation has also demonstrated the physical 
significance of AD which drastically reduces the ac- 
cretion rate compared with the Ohmic-only case. Al- 
though such a comparison may be unfair since our 
Ohmic-only simulation contains net vertical magnetic 
flux, Ohmic-only simulations with similar setup and zero 
net vertical flux generally yield total stress level that is 
only slightly smaller t han ours (e.g., [Turner et all 120071: 
lllgner fc Nelsonll2008f ). Clearly, with AD taken into ac- 
count, zero net vertical magnetic fleld geometry is far 
from being capable of driving rapid accretion in the in- 
ner region of PPDs. 

3.3. Net Vertical Flux Run with Both Ohmic Resistivity 

and AD 

The inclusion of net vertical magnetic flux in our fidu- 
cial run with both Ohmic resistivity and AD (0A-b5) 
makes the field geometry more favorable for the MRI as 
we originally expect. The bottom panel of Figure [2] illus- 
trates the time evolution of Maxwell stress. The initial 
field configuration is unstable to the MRI, which gives 
rise to channel-flow like behaviors and initiates some tur- 
bulent activities in the surface layer. However, we find 
that quite surprisingly, the system then quickly relaxes 
into a non-turbulent state in about 10 orbits with the 
MRI suppressed completely. The laminar configuration 
is then maintained for the remaining of the simulation 
time^. 

As the system settles to a completely laminar state, 
we extract the exact vertical profiles of various physical 
quantities, as shown in Figure [SJ The magnetic field 
is strongest within about ±2.5i7 of the disk midplane, 
where the field is essentially constant due to the large 
resistivity and AD. There is no distinction between active 
layer and dead zone as the entire disk is laminar. The 
disk becomes magnetically dominated beyond z = ±4i7. 
The FUV ionization front is located at about z ~ ztA.bH 
as seen in the sharp increase of Am and Ha profiles, 
which also corresponds to the point where the gas density 
starts to deviate from Gaussian and follow an exponential 
profile. 

Being a 3D time-dependent simulation, the fact that 
the system reaches a steady laminar state suggests the 
stability (particularly, against the MRI) of the config- 
uration. This stability can be qualitatively understood 
using the criterion based on MRI simulations that in- 
clude individual non-ideal MHD effects. In regions near 
the disk midplane where Ohmic resistivity dominates, 
the Ohmic Elsasser number is below one within about 
z = ± 2H, too small for the MRI to operate (|Turner et al.l 
|2007| ). Beyond this region where AD is the dominant 
non- ideal MHD effect, which requires weak magnetic field 
for the MRI to operate, the magnetic field is too strong, 
as judged from Figure 16 of IBai fc Sto^ (|20T1 . Even 
the FUV ionization increases Am substantially beyond 
z = zL4.5H, the disk has already become magnetically 
dominated (/3 < 1) in these regions. The suppression 

^ To justify the validity of this result, particularly that it is not 
due to an unrealistic initial condition, we restart from the end of 
run 0-b5 with AD turned on and find that also in about 10 orbits 
of time, the system settles to the same laminar state. 
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Fig. 5. — Vertical profiles of various quantities in the fiducial run with both Ohmic resistivity and AD and net vertical magnetic flux of 
/3o = 10^ (OA-b5). Upper left: gas pressure and magnetic pressure. Lower left: Elsasser numbers for Ohmic resistivity (A), Hall term (Ha) 
and AD (Am), together with the plasma /3 = Pgas/fmag- Note that wo evaluate Ha while the Hall effect is not included in the simulations. 
Upper right: three components of gas velocity, where the bold green curve is for the vertical velocity. The Alfvcn points are indicated as 
black dots, while open circles mark the base of the outflow. Lower right: three components of the magnetic fleld, where the bold green 
curve is for vertical fleld. 



of the MRI can be understood as a result of magnetic 
field amplification in the surface layer during the initial 
growth from the MRI-unstable configuration. The MRI 
is quenched once the field becomes too strong for it to 
operate because of AD. 

The fact that MRI is suppressed i n th e disk s urface 
layer implies that the hypothesis of iBail ()2011a| ) does 
not always hold, where it was assumed that the mag- 
netic field can be amplified by the MRI to maximize the 
efficiency of the MRI. In other words, we see that the 
magnetic field is amplified to much greater extent that 
the MRI is suppressed. Neverthele ss, this resu lt is not 
inconsistent with the framework of iBail ()2011al) since it 
only estimates the upper limit of MRI-driven accretion 
rate, and complete suppression of the MRI simply means 
the rate is zero. 

The most prominent feature of the laminar state is 
a strong outflow that leaves the vertical boundaries. 



All three components of the velocities exceed the sound 
speed relative to the background Keplerian fiow, with Vz 
reaches about 3cs at the vertical boundaries, and an out- 
flow mass loss rate of pvz = 1.5 x 10~^poCs from each side 
of the simulation box. The Alfven critical point, given 
by Vz = vaz = ^zl \fp^ is contained within our simu- 
lation box and is indicated in the upper right panel of 
Figure El Moreover, according t o the conventional AbI- 
inition ([Wardle fc Koenigllll993() . we define the base of 
the outflow at the location where the azimuthal velocity 
starts to become super-Keplerian, and the it will become 
more obvious later about this definition (see Figure [5]) . 

Our fiducial run 0A-b5 suggests that the structure of 
the disk has a pure one-dimensional (ID) profile. To 
verify this result, we perform a new simulation (run S- 
0A-b5) with everything kept the same except that the 
horizontal domain is reduced to 0.12577 x 0.2577 resolved 
by 4 X 4 cells. We will refer to this type of runs as "quasi- 
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ID" simulations. Such small box is obviously too small 
to study the MRI, but we find that except for different 
initial evolution of the original MRI-unstable configura- 
tion, the system relaxes to exactly the same laminar state 
with a strong outflow as in the fiducial run 0A-b5. More- 
over, we have further checked that although the initial 
evolution involves horizontal variations, such variations 
vanish as the system relaxes to the final steady state^°. 
Therefore, we conclude that the inclusion of AD makes 
the disk structure purely ID. 

In sum, we have seen that the three fiducial simula- 
tions differ with each other in only one piece of physics, 
while show dramatically different characteristics. These 
results best demonstrate the importance of including AD 
in the study of gas dynamics in PPDs, as well as the sig- 
nificant role played by the magnetic field geometry. They 
strongly suggest that the MRI operates extremely ineffec- 
tively in the inner region of PPDs. In the mean time, the 
new laminar solution from our run 0A-b5 strongly points 
to the magnetocentrifugal wind as the most promising 
driving mechanism for PPD accretion. The richness of 
the new findings deserves detailed study, and we devote 
the next section to address the nature of the new laminar 
wind solution. Moreover, the ID nature of the new solu- 
tion allows us to perform the simulations using very small 
horizontal domains which tremendously reduces the com- 
putational cost. 

4. NATURE OF THE WIND SOLUTION 

In this section, we do detailed analysis of our new lam- 
inar wind solution from run 0A-b5. Table [5] summarizes 
the main physical properties of our fiducial laminar wind 
run together with a number of companion runs as they 
will be elaborated in the text. 

4.1. Field Line Geometry and Wind Launching 

We first consider the geometry of poloidal magnetic 
field lines, from which much insight can be gained on 
the wind launching mechanism. Using the magnetic field 
vectors from our run 0A-b5, we integrate a poloidal field 
line from the disk midplane all the way to the vertical 
boundary of our simulation box, and show the results in 
Figure [51 In the mean time, we overplot the direction of 
velocity vectors as red arrows. 

The field line is straight within about ±2.5H from the 
midplane due to the extremely large resistivity, where 
the gas and magnetic field are essentially decoupled as 
we discussed in the previous section. The field lines start 
to bend once the gas become partially coupled to the 
magnetic field, characterized by Ohmic as well as AD El- 
sasser numbers A and Am exceeding unity, which occurs 
a.t z zL2.3H, as can be seen from the bottom left panel 
of Figure \5\ We label this point as the launching point 
in Figure [51 Beyond this point, the azimuthal magnetic 
field decreases rapidly with height, creating large current 
density in the radial direction. Together with the verti- 
cal field, we obtain the following force balance equation 

However, if we start from a pure one-dimensional profile, the 
system will remain one-dimensional but does not relate to a steady 
state. 




x/H 



Fig. 6. — The poloidal field line geometry in our fiducial run 
OA-b5 (blue solid line). Ovcrplotted are the unit vectors of the 
poloidal gas velocity (red arrows). The location of the wind launch- 
ing point, the plasma fi = 1 point, the FUV ionization front and 
the Alfvon point are indicated (black dash-dotted). Also marked 
is the location at the base of the wind (green dashed). 

in the azimuthal direction 

B^-^~lp^y-^o , (15) 

which states that the Lorentz force is balanced by the 
Coriolis force. This explains the increase of radial veloc- 
ity in the upper right panel of Figure [S] beyond about 
zL2.5H, as well as the direction of the arrows in Figure 
[6] at the same locations. In this region, AD is the domi- 
nant non-ideal MHD effect. The radial motion of the gas 
makes the velocity vector deviate from the direction of 
the magnetic field, which drags and bends the magnetic 
field lines towards the same direction via the ion-neutral 
drag. Correspondingly, \Bx\ increases in strength. 

The field lines tend to become straight again as gas 
density decreases and the flow becomes magnetically 
dominated with total plasma /3 < 1. Beyond this point, 
the gas only has limited effect on the fleld line geom- 
etry. Eventually, at the base of the wind (located at 
Zb = ±4.7 H), deflncd as the location where the az- 
imuthal velocity exceeds the Keplerian velocity, the mag- 
netocentrifugal acceleration starts to operate. Through- 
out this paper, we define the region between z = ztz^ as 
the disk zone, while regions beyond as the wind zone. 

We have also indicated the location of the FUV ion- 
ization front, beyond which the gas behaves more or less 
as ideal MHD due to the large ionization fraction. We 
see that the poloidal velocity of the gas is aligned with 
the poloidal magnetic field beyond the FUV front, as ex- 
pected for an ideal MHD wind. Below the FUV front, 
gas velocity deviates from the direction of the magnetic 
fields as a result of AD and Ohmic resistivity. In this 
fiducial run with f3o = W^, the location of the FUV front 
overlaps with the base of the wind, which we note is just 
a coincidence. The importance of the FUV ionization 
will further be discussed in Section [5^ 
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TABLE 2 

Wind Properties from All Fiducial Simulations. 



Run 


"Max 




T-iMax 
















Ep 




Er 




OA-b5 


2.27 X 10" 


-4 


1.07 X 10" 


-4 


3.06 X 10" 


-5 


4.69 


6.10 


1.90 X 10" 


-3 


7.66 X 10" 


-4 


4.38 X 10" 


-4 


S-OA-b5 


2.42 X 10" 


-4 


1.06 X 10" 


-4 


2.86 X 10" 


-5 


4.63 


6.22 


2.02 X 10" 


'3 


8.05 X 10" 


-4 


4.03 X 10" 


-4 


E-OA-b5 


1.38 X 10" 


-4 


1.05 X 10" 


-4 


2.72 X 10" 


-5 


4.56 


6.23 


1.61 X 10" 


-3 


7.30 X 10" 


-4 


4.21 X 10" 


-4 


S-OA-b5-12H 


2.09 X 10" 


-4 


9.67 X 10" 


-5 


4.13 X 10" 


-5 


4.45 


5.33 


1.28 X 10" 


-3 


4.63 X 10" 


-4 


2.28 X 10" 


-4 


S-OA-b5-20H 


2.60 X 10" 


-4 


1.13 X 10" 


-4 


2.23 X 10" 


-5 


4.70 


7.12 


2.77 X 10" 


-3 


1.18 X 10" 


-3 


5.32 X 10" 


-4 


S-OA-b5-24H 


2.58 X 10" 


-4 


1.18 X 10" 


-4 


2.02 X 10" 


-5 


4.83 


7.67 


3.25 X 10" 


-3 


1.14 X 10" 


-3 


8.55 X 10" 


-4 



c«Max: normalized Maxwell stress in the disk zone; T^^": wind stress at wind base; M^j^: total mass loss rate (from both sides of the disk); 
2;,: location of the wind base; Zj\: location of the Alfven point; P^^: rate of work done at the shearing-box boundaries; Ep: energy loss 
rate due to Poynting flux; En: kinetic part of the energy loss rate. All quantities are in natural unit (po = Cs = Q, = 1). 



4.2. Conservation Laws 

A laminar magnetized disk wind is characterized by 
the conservation of mass, angular momentum and energy 
along poloidal streamlines and magnetic field lines (which 
are aligned with each other in ideal MHD). They are very 
useful for diagnosing t he mechanism for wind launching 
and acceleration (e.g., iPelletier fc Pudrit3ll992t ISpruitI 
|1996[ ). In the shearing-box framework, the counterparts 
of the se conservation laws are derived by iLesur et"al] 
((20T1) . It was found that poloidal gas streamlines do not 
necessarily follow exactly the poloidal field lines, which 
introduces new terms in the conservation laws. Never- 
theless, in the absence of strong magnetic field, the devia- 
tions are small (as we checked in our simulations) there- 
fore, we just consider conventional terms in the treat- 
ment. 

Starting from mass conservation, we have 



const kB, 



(16) 



at each side of the disk. In practice, because we con- 
stantly add mass to the simulation box to maintain 
steady state, we find that k does not become constant 
until beyond z ^ ±3H for our run 0A-b5. 

Specific angular momentum is conserved along stream- 
lines, as long as k is constant: 



(17) 



where C 



fix/ 2 is the fluid part of the specific 



angular momentum. The partition between C and By/n 
describes the angular momentum exchange between gas 
and magnetic field. 

Energy conservation is expressed by the Bernoulli in- 
variant in the ideal MHD regime, and we expect it to 
strictly hold only beyond the FUV ionization front (at 
z « ±A.7H). It reads 



i?Bcr — Ek 



Et + EA; 



Ep 



^ + cl log(p) + - 



Byv; 



(18) 



where i^Bcr represents the specific energy along a stream- 
line, with the four terms denoting kinetic energy, en- 
thalpy, potential energy and the work done by the mag- 
netic torque, respectively, and 



Note that the full velocity u (rather than v) enters Ek, 
and the potential energy for the shearing-box is 



(20) 



KBy/p 



(19) 



6 = --rja;2 4. if|^2 
^ 2 2 

To test these conservation laws, we integrate the veloc- 
ity profiles shown in Figure [5] to obtain the streamline. 
We start from z = 2H where k approaches the constant 
value and trace the streamlines to the vertical boundary 
z = 8H, setting x = at z = 2H (note that by construc- 
tion the conservation laws do not depend on the choice 
of the zero point of x). This range covers the most in- 
teresting regions of the disk where wind launching and 
acceleration take place. 

In Figure [7] we show the vertical profiles of individ- 
ual terms in the specific angular momentum and energy 
for streamlines obtained from run 0A-b5. On the left 
panel, we see that the black dashed line flattens beyond 
z « 2.5H, indicating angular momentum conservation. 
This is consistent with our expectations (/ is conserved 
once K approaches constant). It is clear that the increase 
of gas angular momentum (|£|) is compensated by the 
magnetic torque. On the right panel, we see that the 
black dashed line flattens beyond z > AH, indicating 
energy conservation. This is again consistent with our 
expectations (i^Bcr is conserved once the gas obeys ideal 
MHD, i.e., beyond the FUV ionization front). The rise of 
the blue curve in the energy plot indicates flow acceler- 
ation, which is mostly compensated by the reduction of 
potential energy (red). This is the basic picture of mag- 
netocentrifugal acceleration, which is not surprising since 
the bending angle of the field lines relative to disk normal 
in the wind zone is about 45°, well above the threshold 
of 30°. In this frame of reference, the magnetic torque 
plays a minor role, although its effect and the centrifugal 
po tential are intercha ngable by shifting the zero point of 
X (|Bai fc Stonell2013[ ). Thermal pressure gradient only 
plays a diminishing role in the acceleration process. 

4.3. Angular Momentum Transport 

Angular momentum transport can be achieved by ra- 
dial redistribution, due to the Rxj) component of the stress 
tensor as discussed in the beginning of Section [31 as well 
as by vertical extraction, due to the zc/j component of the 
stress tensor 

T^^^ = {pVyVz)\±z, , ^^^^ 

E^"""^ = i-ByBz)\±zb , 
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Fig. 7. — Angular momentum (left) and energy (right) conservation along a streamline in the disk sur face from z = 1H to 8iT for the 
laminar wind solution in our fiducial run OA-b5. For the energy plot, the various terms in Equation l lf 811 are represented by: i?Ber (black 
dashed), Ej<^ (blue solid), Ej- (magenta solid), (green solid) and </> (red solid). For the angular momentum plot, the various terms in 
Equation l lf 71 1 arc represented by: / (black dashed), H (blue solid), —By/n (red solid). 

try considerations, which wih further be discussed in the 
next subsection. We see that if Tr^ and T^^ are com- 
parable with each other, than vertical transport would 
be more efficient than radial transport by a factor of 
about R/H. Since we are considering thin accretion 
disks (where shearing-box approximation applies), it is 
in general much easier for disk wind to transport angu- 
lar momentum compared with turbulent transport. 

We start by examining the radial transport of angular 
momentum. Although radial transport is usually a re- 
sult of turbulence which gives correlated fluctuations of 
velocity and magnetic field in the radial and azimuthal 
direction, ordered velocity/field structure also produce 
radial transport. In Figure [51 we show the vertical pro- 
files of the i?0 components of the Reynolds and Maxwell 
stresses. Clearly, even there is no turbulence, a non- 
zero Maxwell stress given by ordered magnetic fields 
dominates the radial transport in the disk zone. This 
corresponds to the u ndead zone scenario discussed in 
[Turner fc SmmI (|2001 . We also sec that contribution 
from the Reynolds stress in the disk zone is completely 
negligible. Integrating the Maxwell stress across the disk 
zone, we find a « aMax ~ 2.3 x 10"''. 

Meanwhile, Figure [5] best demonstrates the advantage 
for dividing the profile of our the laminar wind solution 
into the disk zone and the wind zone: the Reynolds stress 
is by definition zero at the base of the wind, and increases 
rapidly in the wind zone as a result of radially-outward 
and supcr-Kcplerian motion consequent of the magneto- 
centrifugal acceleration; the Maxwell stress has a promi- 
nent peak at the base of the wind, and decreases in the 
wind zone, as magnetic energy is converted to the kinetic 
energy of the outflowing gas. 

For the vertical transport, we measure the z<f> compo- 
nent of the stress tensor T^^ at the base of the wind. 
Again, our choice of the wind base shows its advan- 
tages: by definition, Reynolds stress is zero, and one 
only needs to consider the Maxwell stress, which we find 




Fig. 8. — The vertical profiles of the R<j> components of the 
Reynolds (red dashed) and Maxwell (blue solid) stresses in our 
fiducial run OA-b5. The disk is divided into the wind zone and 
disk zone, given hy Vy = 0. 

where subscripts corresponds to the location for the 
top and bottom sides of the disk, which we choose to 
be at the base of the wind (the reason will become clear 
shortly). The torque from the outflow can be obtained 
by simply multiplying Tztj, by the radius R (which is un- 
specified in shearing-box). 

The total rate of angular momentum transport, assum- 
ing steady-state accretion in the disk zone, can be ob- 
tained by generalizing equation (fT3|) to include the wind 
torque, which gives 



M 



27r 



dzT] 



Rcj} ■ 



in 



2n , 



(22) 



where we have assumed T- 



zcf)\zt, 



-Tz4,\-zt from symme- 
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is T. 



Ma 



1.07 X lO^Voc^ 



Hereafter, we will re- 



fer to the z4> component of the Maxwell stress at the 
base of the wind as "wind stress". We note that since 
T^'^ = —BzB^ where is constant, from Figure [S] 

we see that T^'^^ varies smoothly within the disk and 
should peak between z = zL2.5H. Although wind stress 
eventually drives accretion, a small fraction of the gas 
flow undergoes "decretion" to bend magnetic field lines 
(between the launching point and the wind base, see Fig- 
ure 15]). This outwardly directed gas flow is properly ac- 
counted for by measuring the wind stress at z = ±Zh 
rather than at z ^ ±2.5iJ. 

Applying the MMSN disk model into Equation ([22]) , we 
can estimate the accretion rate driven by the combined 
effect of radial and vertical transports 



0.82 



10- 



R 



-1/2 
■AU 



4.1 



Pocj 



10 



R 



-3/4 



AU 



(23) 



where Af_8 is the accretion rate measured in 10 Mq 
yr~^. We see that radial transport by the Maxwell stress 
(from ordered magnetic field) in the disk zone can drive 
accretion for about 2 x 10~^Mq yr~^. This number is 
already significantly larger than the case for run OA-nobz 
with zero net vertical magnetic flux, yet still too small 
to account for the typical accretion rate for PPDs. The 
accretion driven by wind transport, on the other hand, 
reaches 4.3 x 10~*Mq yr~^, which is sufficient to account 
for the observed accretion rate for most PPDs. 

We see that our laminar wind solution simultaneously 
solves two problems facing the conventional MRI scenario 
and the conventional wind scenario described in Section 
[TJ First, it efficiently drives disk accretion that matches 
the rates inferred from observations, which can hardly 
be achieved in the MRI scenario. Second, it launches 
the magnetocentrifugal wind with only a very weak net 
vertical magnetic field, in contrast with the conventional 
understanding that requires equipartition field at disk 
midplane. Therefore, the rate of wind transport in our 
scenario is only moderate and does not result in rapid 
depletion of the disk. 

4.4. Symmetry and Strong Current Layer 

The previous subsection demonstrates the success of 
the magnetocentrifugal wind scenario as the the driv- 
ing force of disk accretion. However, one problem was 
ignored in the previous discussion. The laminar wind 
solution we obtained obeys the "odd-z" symmetry: 



B,{z) = B,{~z) 

By{z) - By{-Z) 

B,{z) = B,{~z) 



Vy{z) = -Vy{~z) , (24) 

Vz{z) = -Vz(-Z) . 



As a result, the wind field lines on the top and bottom 
sides of our simulation box bend to opposite radial and 
toroidal directions, as illustrated on the left of the Car- 
toon picture in Figure 13 However, physically, one ex- 
pects the field lines on the top and bottom sides of the 
box to bend toward the same direction, pointing away 
from the central star, as illustrated on the right of Fig- 
ure [9] A particularly disturbing fact is that in the case 





odd-z symmetry 



even-z symmetry 



Fig. 9. — Cartoon illustration of the geometry of wind magnetic 
field lines in shearing-box simulations. The shearing-box approxi- 
mation ignores all radial gradients (exeept the shear), henee does 
not contain information about the location of the central object. 
The natural wind geometry launched from shearing-box has the 
"odd-^" symmetry where the field lines on the top and bottom of 
the box bend toward opposite directions (left). In a physical sit- 
uation, where the location of the central object is fixed, the field 
lines on the top and bottom sides of the box should bend to the 
same direction away from the central object. 

with odd-z symmetry, the "wind" does not transport an- 
gular momentum, since Tztj, at the top and bottom sides 
of the disk have the same sign and value, and cancel 
out. The wind angular momentum transport mechanism 
works only for a physical wind geometry. 

We note that all radial gradients except the azimuthal 
shear are neglected in the shearing-box approximation, 
hence one can not tell whether the location of the cen- 
tral star is in the "inner" or "outer" (radial) side of the 
box. In other words, the two radial sides of the box are 
symmetric. In our simulations, we find that the direction 
that the wind field lines bend as we start from the initial 
configuration is totally random, and the top and bottom 
sides of the disk evolve independently^-'^. This suggests 
that the chances to find a solution with the unphysical 
odd-z symmetry and with a physical geometry are equal. 
In fact, our ID run S-0A-b5 results in a solution with 
the physical geometry. 

In Figure [TOl we show the profiles of the magnetic field 
and velocity field from our run S-OA-b5-F3 that have 
the physical wind geometry. We see that this physical 
solution is NOT a solution under the "even-z" symmetry: 



S,(z) = -B.,{~z) 

B.y{z) = -By{-Z) 

Bz{z)^B,{-z) , 



Vx{z) = Vx{~z) , 

Vy{z)^Vy{-z) , (25) 

Vz[z) = -Vz{-Z) . 



which is almost exclusively use d for cons tructing semi- 
analytical local wind solutions (jWardle &: Koenigl 199^; 

We have performed more than one copies of each 3D and the 
quasi-lD fiducial simulations that verifies the randomness. 
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Qgilyie fc Li7ioll2001l: iKonigl at al]|2010l : ISalmeron et"al] 
201 it ). Very interestingly, the physical wind solution fol- 
lows exactly the odd-z symmetry solution we obtained 
before, as can also be seen from Table [2] that all physi- 
cal quantities measured in this run agree with those in 
our 3D fiducial run 0A-b5 within 5%. In particular, the 
wind stress T^'^^ is almost exactly the same. The only 
difference is that the horizontal field lines (and horizontal 
velocities) are flipped at one side of the box to achieve 
the physical wind geometry. This flipping is mediated 
by a sharp transition of the horizontal fields, which ex- 
hibits as a strong current layer. The strong current layer 
is not located at the midplanc, but is located at about 
z = +3H (or at z = —3H, and the selection is random). 
Therefore, the physical solution does not have any sym- 
metry about the disk midplanc. 

The location of the strong current layer roughly corre- 
sponds to where both the Ohmic and AD Elsasser num- 
bers A and Am become greater than 1, as seen from the 
bottom left panel of Figure [5l The reason for such large 
offset from the disk midplane is that magnetic diffusion 
near the disk midplane is so strong that the gas and 
magnetic field are essentially decoupled. Correspond- 
ingly, the magnetic field lines must be straight and cur- 
rent is excluded, and the strong current layer can exist 
only when magnetic diffusion becomes weaker. 

The right panel of Figure [TU] demonstrates the velocity 
profile of the physical wind solution. The strong cur- 
rent layer exhibits as a seemingly small feature at about 
z = +3H in this plot, which is in linear scale. As we 
zoom into this region shown in the inset, we find that 
the gas has a large inflow velocity of about O.Scg at the 
location of the strong current layer. The large inflow is 
directly the consequence of the wind stress: the Maxwell 
stress exerted at the base of the wind is released at the 
strong current layer, driving a large inflow which is es- 
sentially how accretion proceeds. More specifically, it is 
the balance between the Lorentz force and Coriolis force 
that leads to the large inflow and accretion in the strong 
current layer (see equation ([T5|) ). 

To further demonstrate the effectiveness of accretion 
in the current layer, we note that for a MMSN disk, in 
order to have accretion to approach 10~^ Mq yr^^, a 
bulk radial drift velocity of about 10~^Cs is required 

A^_8«2.5(^)i?X^/^ (26) 

In the physical wind solution, we find that only a tiny 
fraction (about ^ 5 x 10~^) of disk mass is contained 
in the strong current layer, but it drifts at very large 
velocities (^ 0.3cs). The combined effect is an efficient 
accretion that easily accounts for the typical accretion 
rates observed in PPDs. More specifically, by integrating 
the radial mass flux (/ pvxdz) across the disk zone, we 
find the mean radial velocity of the gas to be —1.69 x 
10~^Cs. Plugging to Equation ([25)1 . this would lead to an 
accretion rate of 4.2 x 10~^A/q yr~^, which agrees very 
well with the estimate in Section 

Finally, we note that the strong current layer is not a 
thin current sheet, since at its location there is still sub- 
stantial magnetic diffusion (mainly AD). It has a thick- 
ness of about 0.25H, which is well resolved in about 6 



cells in our quasi-lD simulations. We have found that 
the strong current layer is stable in our quasi- ID fiducial 
run S-0A-b5. The stability of the strong current layer 
in 3D requires justification in the more geometrically ap- 
propriate global simulations^^. 

4.4.1. Wind Solution from Even-z Symmetry 

We have discussed that the odd-z symmetry is the most 
natural symmetry in shearing-box but it is unphysical for 
a disk wind. A physical wind solution can be achieved by 
flipping the horizontal field and stream lines at one side 
of the disk, which yields a solution that contains a strong 
current layer offset from the disk midplane. Because the 
two solutions have exactly the same properties except the 
flip, it is natural to classify them as the same solution. In 
this subsection, we further consider a solution that obeys 
even-z symmetry. 

We design a new 3D simulation E-0A-b5 which is ex- 
actly the same as our fiducial run 0A-b5 except that we 
use only half of the shearing-box, with the lower bound- 
ary located at the disk midplane. At the lower bound- 
ary we use the standard reflection/conduction boundary 
condition so that even-z symmetry is enforced. We find 
that after a brief period of the MRI, the system also set- 
tles into a pure laminar state, with velocity and magnetic 
field proflles shown in Figure[TT] (where we have recovered 
the other side of the disk based on even-z symmetry). 
Under even-z symmetry, the horizontal magnetic field at 
the midplane is enforced to be zero. Due to the large 
resistivity (which can not support a significant current), 
the field is kept very close to zero within z ~ ±27J where 
both Ohmic and AD Elsasser numbers are much less than 
1. Beyond z = zt2H, the horizontal field increases, and 
then beyond z ~ ±3i?, the fleld configuration recovers 
to almost exactly the same solution we obtained under 
odd-z symmetry. 

We note that the fleld configuration similar to our 
even- z syrn metr y solution has previously been reported 
byO (jl996| ) and (Wardlel()1997t ). and it was discussed that 
in the presence of low conductivity in the midplane, wind 
solution can be constructed with midplane magnetic fleld 
pressure much smaller than the gas thermal pressure. 
Our solution conflrms their findings and demonstrates 
it robustness and uniqueness since it is constructed in a 
self-consistent and evolutionary manner. 

The velocity profiles from our even-z symmetry solu- 
tion are exactly the same as the odd-z symmetry solu- 
tion in the wind zone (except the flip). In the disk zone, 
we see that the radial velocity has two dips located at 
z ±3H. The dips are much shallower (only about 
O.OlCs) than the dip due to the strong current layer in 
our run S-0A-b5, but they are much wider. Integrating 
the pVx over height, we obtain the radial mass flux to be 

^■^ We have repeated our 3D fiducial run OA-b5 with a different 
random seed which evolves to a physical wind solution with a strong 
current layer. We find that the strong current layer is maintained 
for about 100 orbits but eventually escapes the simulation box 
and the odd-z symmetry is recovered. This is likely to be due 
to the intrinsic limitations of the shearing-box, where curvature 
terms are neglected and hence favors the odd-z symmetry solution. 
Moreover, as the fast magnotosonic point is beyond our simulation 
box, the stability of the strong current layer can also be affected 
from outside of the simulation domain. 
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Fig. 11. — Same as Figure [TOl but for run E-OA-b5 where an even-2 symmetry ||25| I is enforced (see Section |4.4.1| I. The inset on the 
right panel shows the zoomed-in view of the radial velocity profile in the disk zone. 



— 1.67 X 10~^Cs, which is almost exactly the same as in 
run S-0A-b5. This is not surprising since the properties 
of the wind in the two solutions are exactly the same. 

In sum, wc find that the suppression of the MRl and 
wind launching is inevitable in the inner region of PPDs 
around lAU, and the property of the solution in the wind 
zone is independent of the large-scale field geometry (or 
symmetry) . 

4.5. Outflow 

Strong outfiows are launched in our fiducial simula- 
tions, and we have measured that the rate of mass out- 
flow from each side of the box is about 1.4 x 10~^poCs 
from the quasi-lD run and 1.5 x 10~^poCs from the 3D 
run. In Table [21 we also list total mass loss rate, 

Mn, = \pVz\hot + \pVz\top (27) 

as the sum of the mass loss rate from the two vertical 
boundaries for our fiducial runs. Although mass loss is 
not significant for the duration of our simulations, the 
value is relatively high that without replenishing to the 



disk, the mass loss timescalc is only 8 x lO^il"^, which 
amounts to only about 10^ years. Moreover, the mea- 
sured wind mass loss rate is in fact comparable to the 
mass accretion rate driven by the wind itself, which is 
inconsistent with the observationally inferred ratio that 
the mass loss rate i s only about 10% of th e accretion rate 
(jCabrit et al.lll990l: iHartigan et al.lll995[ ). 

We note that, however, the magnetocentrifugal wind is 
intrinsically a global phenomenon. A global wind solu- 
tion is not fully determined until all critical points (or 
surfaces), namely, the slow magnetosonic, Alfven and 
fast magnetosonic points are passed, where the fast mag- 
netosonic point is far beyond the extent of our simula- 
tion box. This leads to one extra degree of freedom in 
the wind solution. Moreover, the global structure of the 
wind, and the location of the critical points, are set by 
the interplay between the microphysics in the disk, and 
the large-scale field structure (i.e., the magnetic flux dis- 
tribution over the entire disk), where the latter is not 
reflected in local shearing-box simulations. Therefore, 
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while useful for studying wind-launching, the shearing- 
box approximation has its limitations and the measured 
rate of outflow from our simulations should only be taken 
as reference and it may significantly overestimate the 
mass outflow rate in reality. 

To further clarify the discussion above, we perform 
three additional quasi-lD simulations, changing the box 
size to 12H, 20H and 2411 respectively, and the three 
runs are labeled as S-OA-b5-12H, S-OA-b5-20H and S- 
OA-b5-24H. We measure a series of quantities discussed 
before and also list the results in Table [2l We find that 
as we increase the height of the simulation box, the lo- 
cation of the Alfven point moves to higher altitudes. In 
the mean time, the mass loss rate decreases. In fact, 
these two quantities are connected, and in shearing-box, 
the d efinition of the Alfven point indicates ()Bai &: Stond 
[20T3h 



(28) 



where \a denotes the value at the Alfven point. There- 
fore, higher location of the Alfven point indicates lower 
density, hence smaller mass loss rate. 

The trend of deceasing mass outflow rate with vertical 
box size reflects one limitation to shearing-box: the grav- 
itational potential increases quadratically with vertical 
height, which increasingly hinders the gas from escap- 
ing. The outflow in our simulations is possible because 
gravitational potential is cut off at the vertical bound- 
aries. This limitation of shearing-box has been noticed 
and discussed in a num ber of recent works on ideal MHD 
simul ations of the MRI (jFromang et al.l2012l : iLesur et al.l 
l2012f ). Because the depth of the shearing-box potential 
is about the same as t he re al depth of the potential at 
z ~ ±i?, iBai fc Stond ()2013[ ) suggest that the mass loss 
rate can be roughly estimated by performing shearing- 
box simulations with a box height of Lr ~ 2R, or to 
study the height dependence of M^, and extrapolate the 
resuh to - 2R. In MMSN, we have H/R - 0.03 at 
1 AU (note that H/R is not specified in the shearing- 
sheet approximation). Such a thin disk would require 
a shearing-box height of ^ 65H, which is numeri- 
cally unpractical. Based on our three quasi-lD runs with 
Lz — 12H to 24iJ, the trend can be roughly described 
by ''^ l/iz- Extrapolating this trend, we obtain a 
rough estimate of mass loss rate to be ^ 6 x 10~^ 
in code unit. This level of brings the mass loss time 
scale to about 10^ years, which is about an order of mag- 
nitude longer than our original naive estimate. Although 
still somewhat short, the mass reservoir in the outer disk 
(which contains much more mass than that in the inner 
disk) may be sufficient to feed the inner disk to compen- 
sate for the mass loss. 

It is also interesting to note that although the mass 
outflow rate is largely affected by the vertical box size, 
the angular momentum transport depends on the vertical 
box size very weakly. For box sizes of 12H, IQH, 20H 
and 24iJ, the values of T^^'^ at the base of wind are 

found to be 0.97 x lO"'', 1.06 x lO""*, 1.13 x lO""* and 
1.18 X 10"'' respectively, hence increasing box size leads 
to more or less the same (but slightly stronger) wind- 



driven accretion. 

In sum, the trend that the mass outflow rate decreases 
with box size and mass accretion rate increases with box 
size implies that in a real system, the ratio of mass ac- 
cretion rate to mass loss rate would be a factor of several 
larger than that obtained in our simulations, and would 
be more consistent with observations. 

4.6. Energetics 

Although our simulations assume an isothermal equa- 
tion of state where total energy is not conserved, it is 
important to verify that our simulation results are ener- 
getically feasible and to examine the energy balance for 
real situations. 

In shearing-box simulations, the energy is injected due 
to the work done by the shearing-box bou ndaries, the 
rate of which is given by (IBai fc Stond[20T3h 

1 /■ 3 3 
-Psh = -j-j^ / -^{pVxVy - BxBy)nL^dydz = -Sc^^Ja , 

(29) 

where the integral over z is performed within z = ±Zf, 
(the disk zone) . Energy is lost through the open vertical 
boundaries, given by 



E = Ek + Ep ^ pl— + ci\v-{vxB)xB n 

. J J z=±zb 

(30) 

where n is the unit vector pointing away from the disk 
in the vertical direction, and we sum over contributions 
from the top and bottom of the wind base. The flrst 
terms in the bracket represent the kinetic energy loss 
and the PdV work done by the mass outflow, and the 
second term represents energy loss from the Poynting 
flux. We have ignored the energy loss term associated 
with internal energy loss due to the mass outflow, since 
we keep feeding mass to the system which balances this 
term exactly. 

In Table [2] we also list the values of Psh, Ex and Ep^ 
normalized in natural units. We see that the sum of E^ 
and Ep comprises of about 60% of the total work done by 
the shear Pgh; hence energy conservation is not violated. 
The extra energy is likely to be radiated away in real 
systems. 

4.7. Radial Drift of Magnetic Flux 

Globally, poloidal magnetic flux drifts radially at ve- 
locity in the presence of a net toroidal electromotive 
force (EMF, E) 

^ ^ . (31) 

Negative vb would lead to accumulation of magnetic flux 
in the inner disk while positive vp would make the disk 
lose magnetic flux to large scales. Moreover, 8y should 
be constant with height so that magnetic flux drifts uni- 
formly across the disk. In global models, the radial drift 
velocity vp is determined b y the radial distribution of 
magnetic flux (|Teitleiir2011[ ). In local models, it is often 
chosen as a fr ee parameter in local she aring-box models 
of disk winds (jWardle fc Koeni^ll993f ) due to the extra 



18 



X.-N. Bai & J. M. Stone 



0.2 




i> 




' 1 




0.15 




1 ' 
, 1 
> 1 




1 

1 - 




0.1 




> 1 
1 1 
1 1 




I - 
1 - 


- 


0.05 


\ 


1 1 




1 
1 




C 

LU 




1 V 
J ^^ 




i. 










V 1 


- ^ ^ ' 
N 


-0.05 




' ; 




1 1 
1 1 


\ 


-0.1 




! 1 
' 1 
' 1 




1 1 

1 1 
1 , 




-0.15 




I 1 
1 1 




1 , 

1 1 








ii 




1, 

1 \ 




i -6 


-4 -2 





2 4 


6 8 



TABLE 3 

Summary of All Simulations in the Parameter Study. 



z/H 

Fig. 12. — The toroidal EMF profile in our fiducial run OA- 
b5. Blue dashed: inductive EMF (f^^); red dash-dotted: AD EMF 

(f^); black solid: the total EMF {By). The Ohmic EMF is very 
close to and is not plotted. The EMFs are normalized to B, and 
are in natural unit of the sound speed Cs. 

degree of freedom discussed in Section 14.51 Our sim- 
ulations in principle have one extra degree of freedom 
because the fast magnetosonic point is not contained in 
the simulation box. We expect that the value of vb in 
our simulations has direct correspondence with this extra 
degree of freedom. 

The toroidal component of the EMF in our simulations 
is the sum of the inductive EMF (f^), Ohmic EMF (£:°) 
and AD EMF which reads 



(32) 



In Figure [12] we show the radial profiles of the three 
toroidal EMFs and their sum, all normalized to B^. We 
see that in most regions of the disk, the inductive EMF 
is essentially exactly balanced by the AD EMF. We did 
not plot the Ohmic EMF because it is always negligibly 
small due to the very small current near the midplane 
and the very small resistivity in the upper layer. The 
sum of all EMFs, is consistent with zero at all locations 
of the disk. The little sharp peaks of the total EMF at 
around z « ±4.7i? are purely numerical features due to 
the sharp variations of the AD coefficient and current 
density at the FUV ionization front. 

Therefore, the wind solution we obtained from 
shearing-box simulations is a zero radial drift solution 
Vb = 0, which corresponds to a stationary magnetic con- 
figuration in the global picture. A direct consequence of 
cb = is that poloidal field lines and poloidal stream- 
lines follow each other in the absence of non-ideal MHD 
terms, as clearly shown in Figure El The fact that vb = 
in our simulations is likely to be a result of our vertical 
boundary condition, where there is no vertical gradient of 
magnetic fields. We note that in the ideal M HD shearing- 
box si mulations of wind launching process bv lLesur et al.l 
()2012|) . significant deviation between poloidal field lines 
and streamlines are found (see their Figure 3). This is 
likely to be a result of a different vertical boundary condi- 
tion: poloidal field is forced to be vertical at ghost zones, 
resulting in a current sheet at the vertical boundary. 



Run 


S/Smmsn 


/3o 


Spuv 




S-OA-b3 


1 


10^ 


0.03 


10"* 


S-OA-b4 


1 


lO* 


0.03 


10-* 


S-OA-b6 


1 


106 


0.03 


10-* 


S-OA-Fl 


1 


10^ 


0.01 


10-* 


S-OA-FIO 


1 


10^ 


0.1 


10-* 


S-OA-nogr 


1 


10^ 


0.03 





S-OA-grOl 


1 


10^ 


0.03 


10-=* 


S-OA-M3 


3 


3 X 10^ 


0.1 


10-* 


S-OA-M03 


0.3 


3 X 10* 


0.1 


10-* 



5. PARAMETER STUDY 

To assess the importance of various physical effects on 
the properties of the laminar wind solution, we conduct a 
parameter study by surveying four parameters and com- 
pare the results from our fiducial model. All simulations 
in this parameter study are quasi-lD. First, we vary the 
vertical net magnetic flux, and consider /3o = 10^, 10^ 
and 10^. These runs are labeled as S-OA-bn, where bn 
denotes /3 — 10". Second, we vary the penetration depth 
of the FUV ionization, and consider Spuv = 0.01 g cm-^ 
and 0.1 g cm-^, and the run labels are attached with 
"Fl" and "FIO" respectively. Third, we consider the ef- 
fect of grain abundance, and consider the case with grain- 
free chemistry (whose run label is attached with "nogr" ) 
and the case with grain abundance egr = 10"'^ (run la- 
bel attached with "grOl" , indicating depletion factor of 
0.1), 10 times the value in our fiducial run. Finally, we 
vary the surface density of our disk model, and consider 
S = 3SMMSN and S = 0.3Smmsn- In the mean time, we 
maintain the same field strength in physical unit, thus 
/3o in these two runs are set to 3 x 10^ and 3 x 10'' re- 
spectively. These two runs are attached with "M3" and 
"M03" , respectively. A list of all these runs are given in 
Table [1 

All the additional quasi-lD simulations are run to 
t ~ 120017"'. We extract the vertical profiles by time 
averaging the second half of the runs and perform the 
same analysis as we did for the fiducial run. In partic- 
ular, we identify the locations of the wind base and the 
Alfven point, the a parameter of the Maxwell stress in 
the disk zone, the zcj) component of the Maxwell stress at 
wind base, and so on. The results of these measurements 
are listed in Table SI Before discussing these physical ef- 
fects in details, we plot in Figure [13] the rate of mass 
outfiow and the wind stress as a function of net vertical 
field strength. 

5.1. Effect of Net Vertical Magnetic Flux 

The vertical net magnetic field is a crucial ingredient 
to launch a magnetic outflow. It allows the field lines 
to be connected to infinity, paving way for the outflow 
to escape from the disk. We see from Figure [T3] that 
the rate of mass outflow increases rapidly with in- 
creasing net vertical field. A fit to the scattered plot 
gives M^u oc Pq^-^"^. In other words, we roughly have 
(X Bz- This result is con sistent with the ideal MHD 
shearing-box simulations by ISuzuki fc Inutsukal (|2009l ) , 
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TABLE 4 

Wind Properties from All Simulations in the Parameter Study. 



Run 








T^Max 




M 








^A 


^sh 




Ep 




Ek 




S-OA-b3 


1.19 X 


10" 


2 


3.37 X 


10- 


3 


4.76 X 


10- 


4 


3.95 


>8.00 


6.93 X 


10- 


2 


4.29 X 


10" 


-3 


3.52 X 


10- 


-2 


S-OA-b4 


1.27 X 


10" 


3 


5.90 X 


10- 


4 


8.70 X 


10- 


5 


3.92 


7.47 


9.31 X 


10- 


3 


2.42 X 


10" 


-3 


3.34 X 


10- 


-3 


S-OA-b5 


2.42 X 


10" 


4 


1.06 X 


10- 


4 


2.86 X 


10- 


5 


4.63 


6.22 


2.02 X 


10- 


3 


8.05 X 


10" 


-4 


4.03 X 


10- 


-4 


S-OA-b6 


3.01 X 


10- 


5 


2.86 X 


10- 


5 


1.00 X 


10- 


5 


4.52 


5.52 


4.62 X 


10- 


4 


2.14 X 


10- 


-4 


4.17 X 


10- 


"5 


S-OA-Fl 


1.91 X 


10- 


4 


9.16 X 


10- 


5 


1.62 X 


10- 


5 


4.17 


7.11 


1.61 X 


10- 


3 


6.14 X 


10- 


-4 


2.81 X 


10- 


-4 


S-OA-FIO 


1.62 X 


10- 


4 


1.66 X 


10- 


4 


5.26 X 


10- 


5 


4.20 


5.52 


2.60 X 


10- 


3 


1.09 X 


10- 


-3 


5.39 X 


10- 


-4 


S-OA-grOl 


1.32 X 


10- 


4 


1.04 X 


10- 


4 


2.69 X 


10- 


5 


4.55 


6.30 


1.63 X 


10- 


3 


7.83 X 


10- 


-4 


3.92 X 


10- 


-4 


S-OA-nogr 


4.96 X 


10- 


4 


1.16 X 


10- 


4 


3.84 X 


10- 


5 


4.92 


5.86 


2.87 X 


10- 


3 


9.03 X 


10" 


-4 


4.62 X 


10- 


-4 


S-OA-M3 


8.04 X 


10- 


5 


3.64 X 


10- 


5 


1.06 X 


10- 


5 


4.89 


6.23 


6.50 X 


10- 


4 


2.67 X 


10- 


-4 


1.38 X 


10- 


-4 


S-OA-M03 


7.34 X 


10- 


4 


3.51 X 


10- 


4 


9.17 X 


10- 


5 


4.38 


6.09 


6.60 X 


10- 


3 


2.79 X 


10- 


-3 


1.36 X 


10- 


-3 



Note: we have duplicated run S-OA-b5 as a standard reference. 
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Fig. 13. — The rate of mass outflow (left) and the zcfi component of the Maxwell stress at the base of the wind (right) as a function 
of vertical net field strength, indicated as l//3o for all our simulation runs. Particularly, blue squares represent simulations with varying 
vertical field strength, red diamonds represent simulations with varying FUV penetration depth, and green circles represent simulations 
with varying grain abundance. 

find the radial transport remains to play a minor role to 
the total accretion rate for all value of /3o explored so far, 
but its contribution increases with increasing net flux. 

Finally, we find that for relatively strong net vertical 
field with /?o ^ 10^, it is less likely to obtain a laminar 
solution with the physical wind geometry: the strong 
current layer would escape from our simulation box and 
one recovers the undesired solution with odd-z symme- 
try. This observation may suggest that maintaining the 
stability of the strong current layer becomes more diffi- 
cult for larger field strength, although the issue can only 
be fully resolved by performing global simulations. 



as well as more recently by iBai fc Stond ()2013l) . Al- 
though the wind launching is associated with vigorous 
MRI turbulence in ideal MHD, the scaling of the mass 
loss rate with vertical background magnetic field is very 
similar. 

On the right panel of Figure fT3l we see that the rate of 
angular momentum transport, characterized by the wind 
stress, increases even more rapidly with background ver- 
tical field strength. A simple fit returns T^^'^ cx /^q^"'^. 

Since accretion rate Ma oc T^'^^, increasing vertical net 

field reduces the ratio of M^/Ma- This is accompanied 
by the outward shift of the location of the Alfvcn point 
with increasing net vertical flux, and for /3o = lO'^, the 
Alfven point is beyond our simulation box. We also note 
for relatively large net vertical flux of /3o — 10^, the wind- 
driven accretion rate (Equation would become too 
large compared with observations. Therefore, a reason- 
able wind-driven scenario in PPDs should involve only 
weak vertical background field with /Jq of the order 10^ 
to a few times 10^. 

Increasing the net vertical field also leads to rapid in- 
crease in the midplane magnetic field (undead zone), 
which leads to larger aMax- Applying Equation (|23p . we 



5.2. Effect of FUV Ionization 

The FUV ionization is another crucial ingredient of 
wind launching: the large ionization fraction in the FUV 
layer makes the gas and magnetic field be strongly cou- 
pled to each other so that it is essentially in the ideal 
MHD regime. The strong coupling between the gas and 
magnetic field is essential for effectively loading mass 
onto open magnetic field lines for magneto-centrifugal 
acceleration to take place. Therefore, we expect the rate 
of mass outfiow to strongly depend on the penetration 
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depth of the FUV ionization. Indeed, we see from Fig- 
ure [13] that increasing Spuv by a factor of 10 results 
in a factor of more than 3 increase of the mass outflow 
rate. Meanwhile, increasing Sfuv also leads to a moder- 
ate increase of the wind stress, as seen on the right panel 
of Figure [131 Correspondingly, the ratio of M^jMa in- 
creases with Efuv> in parallel with the lowering of the 
the Alfven point. 

5.3. Effect of Grain Abundance 

In th e MRI-driven accretion scenario studied in [Bail 
()2011a|) , the predicted accretion rate sensitively depends 
on the size and abundance of grains. In the wind-driven 
accretion scenario, we find that the dependence on grains 
is much weaker. For the outflow rate, there is only a 50% 
difference between the grain-free chemistry and the case 
with 0.1/zm grains with the abundance of egr = lO"-^. 
The wind stress, on the other hand, is almost indepen- 
dent of grain abundance. The weak dependence on the 
grain abundance is mainly due to the fact that the wind 
launching process mainly depends on the gas conduc- 
tivity in the surface layer of the disk where grains only 
play a minor role: either the FUV ionization dominates 
(which is independent of grain abundance), or the long 
recombination time in the low gas density leads to large 
ionization fraction that well exceeds the grain abundance 
(hence grains have only ve ry limited e ffect on the ioniza- 
tion level, see Figure 1 of lBal[20TIal) . The grain abun- 
dance does have a strong effect on the radial transport, 
as we see from Table [2] that ofMax increases relatively 
rapidly with decreasing grain abundance. Nevertheless, 
the overall picture is unchanged since radial transport 
only plays a minor role in driving disk accretion. 

5.4. Effect of Disk Mass 

The surface density of the inner PPDs is highly uncer- 
tain and have not been well constrained observationally 
(situation is better for the outer disk with millimete r 
interferometric observations, e.g.. [Andrews et "aTl 120091 ). 
The MMSN surface density serves as an initial guess and 
is likely on the lower end, while m uch more m assive in- 
ner disks are also speculated (e.g.. [Desch[[2007l) . On the 
other hand, the surface density should decrease with time 
as the disks evolve. Therefore, we further consider disk 
models that are three time (with label "M3") and 0.3 
times (with label "M03" ) the mass of the fiducial MMSN 
disk, while we have also varied the net vertical flux in 
code unit so that the physical magnetic fluxes are the 
same for all runs. Because we have adopted midplane 
density po = 1 as our natural unit whereas in reality 
Po oc S, hence some conversion is necessary. In Table [31 
we see that the wind fluxes T^^'^ in the two cases vary by 
a factor of about 10, and the fiducial case lies in between. 
This means that in physical unit, the stress is more or 
less independent of the disk surface density. The wind 
mass loss rate M^, behaves similarly and is roughly in- 
dependent of S. Therefore, the accretion and mass loss 
rate in PPDs is largely determined by the poloidal mag- 
netic flux distribution through the disk, regardless of disk 
mass. 

6. DISCUSSIONS AND CONCLUSION 



6.1. Summary 

In this paper, we have performed three-dimensional 
shearing-box simulations to study the gas dynamics in 
the inner region of a PPD (at 1 AU). Non- ideal MHD 
effects, namely the Ohmic resistivity, Hall effect and am- 
bipolar diffusion (AD) are crucial in the weakly ionized 
gas as in PPDs. They affect the coupling between the 
gas and magnetic field in different ways and dramati- 
cally change the MHD stability of the system. Partic- 
ularly, Ohmic resistivity dominates in the midplane re- 
gion while AD dominates in the disk surface layer, and 
the Hall dominated regime lies in between. For the first 
time, we have included both Ohmic resistivity and AD 
in vertically stratified simulations. The diffusion coeffi- 
cients are obtained self-consistently by interpolating from 
a pre-computed look-up table (as a function of density 
and ionization rate at fixed temperature) based on equi- 
librium chemistry. 

We first showed that the conventional picture of lay- 
ered accretion fails when AD is taken into account. With- 
out considering AD (the conventional scenario), the MRI 
drives vigorous turbulence in the disk surface (active) 
layer that easily accounts for the observed accretion rates 
in PPDs. However, we find that the MRI activity is 
severely suppressed by AD, and the gas dynamics also 
sensitively depends on magnetic field geometry. More 
specifically, in the zero net vertical flux field configura- 
tion, the disk exhibits extremely weak turbulent activ- 
ity, with total stress parameter a ~ 3 x 10~^, mostly 
contributed from a thin layer above the FUV ionization 
front. The resulting steady state accretion rate is about 
three orders of magnitude too small compared with ob- 
servations. The main reason for such inefficiency is that 
MRI in the AD dominate d regime disfavors z ero net ver- 
tical fiux field geometry (|Bai &: Stond l201ll ). By con- 
trast, in the presence of net vertical magnetic flux, the 
MRI turbulent activity disappears and the flow is com- 
pletely laminar. In the mean time, the disk launches a 
strong wind, which is accelerated to super- Alfvcnic ve- 
locities within our simulation domain. Instead of MRI 
driven accretion in the conventional scenario, angular 
momentum is carried away by the magnetocentrifugal 
wind. 

We provide detailed analysis of the laminar wind so- 
lution. The wind launching process is mainly assisted 
by the vertical gradient of toroidal magnetic field, start- 
ing from about z ^ zt2H where gas and field lines are 
marginally coupled to each other (near the midplane they 
are essentially decoupled). At the base of the wind, the 
azimuthal gas velocity exceeds Keplerian, and the mag- 
netocentrifugal mechanism takes place and efficiently ac- 
celerates the wind flow. We flnd that the magnetocen- 
trifugal wind very effectively carries away angular mo- 
mentum from the disk. Even with a small net vertical 
magnetic flux of /3o = 10^, the launched wind is already 
sufficient to account for the observed accretion rate of 
about 10~^Mq yr~^. The rate of the mass outflow 
is not a well characterized quantity in shearing-box simu- 
lations, which decreases with increasing vertical domain 
size. We speculate that realistic mass loss rate is much 
smaller than that measured in our simulations based on 
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studying the dependence of Af^, on box height; on the 
other hand, the wind stress, hence the wind-driven ac- 
cretion rate increases only shghtly with increasing box 
height, hence our measurement is more rehable. 

The natural symmetry of the wind solution obtained 
in our shearing-box simulations has an odd-z symme- 
try which is inconsistent with a physical wind geometry. 
Meanwhile, we find that the system has an equal prob- 
ability to evolve into an exactly same solution but with 
horizontal velocity and magnetic field flipped at one side 
of the box. The flip takes place in a strong current layer 
with a thickness of ^ 0.25H that is offset from the disk 
midplane (due to the large resistivity at midplane). This 
second solution has the desired physical geometry of the 
magnctoccntrifugal wind, and in this solution, the entire 
accretion flow is carried within the strong current layer, 
with radial velocity achieving ~ 0.1 — 0.3cs, while the rest 
of the disk remains static. We further show that a so- 
lution with enforced even-z symmetry gives exactly the 
same magnctoccntrifugal wind in the disk upper layer, 
which further vcrifles the robustness of wind-driven ac- 
cretion scenario. 

We have performed a parameter study where we show 
that the rate of wind angular momentum transport T^'^^ 
increases rapidly with increasing net vertical flux, with 
approximately T^'''' cx /3o"°-^ Similarly, cx /3o"° ^ 
Moreover, both the wind stress and the wind mass loss 
rate depend only on the physical vertical field strength 
and very weakly depends on the surface density of the 
disk. The FUV ionization which makes the gas in 
the wind zone behave as ideal MHD, is essential for 
wind mass loading. The depth of the FUV ionization, 
which depends on the extent of dust attenuation and 
self-shielding, also affects the efficiency of wind trans- 
port of angular momentum and wind mass loss rela- 
tively strongly. Finally, we find that the grain abun- 
dance, which has been found to play a crucial role in 
determining the extent of the acti ve layer in th e conven- 
tional layered accretio n scenario (jllgner fc Nelson 200^ 
iBai &: Goodman! |2009| ) . has almost no influence on the 
wind scenario. This is mainly because that the wind 
launching process takes place in the upper layer of the 
disk where the chemistry is less affected by the grains 
(since ionization fraction well exceeds grain abundance). 

6.2. Implications 

The fact that the inner region of PPDs is likely to 
be laminar with purely wind-driven accretion processing 
through a strong current layer offset from the disk mid- 
plane may have profound implications on many aspects 
of planet formation, particularly on the following two as- 
pects. 

The laminar inner disk is likely to become the mostly 
favored spot for grain growth and settling, as well as 
planetesimal formation and growth. In the conventional 
picture of layered accretion, the Ohmic dead zone re- 
gion is found to be not completely static, but has ran- 
dom motion due to the sound waves launched from the 
active layers fFleming & Stonc"2003"; 'Oishi fc M ac Low! 
[2009; Turner ct al. 2010; Okuzumi & Hirosc 2011). The 
level of the random motion, albeit much smaller than 



that in the conventional active layers, can still be large 
enough to prevent the solids from settling completely. 
The random gas velocities plays an important role on 
the collision velocit y between dust grains that would im- 
pede grain growth ( Ormel fc Cuzzil[2007[ ). Moreover, the 



stochastic gravitational torque from density fluctuations 
would excite the random velocities among existing plan- 
etesimals to sufficiently large velocities that their mu- 
tual collision s may lead to f ragmentation rat her than 
net grow th dlda et all l2008t [Yang et alj |2012[). More 
recent ly, iGressel et al.l ( 2012 ) and lOkuzumi fc Hirosd 
(|2012t ) found that increasing the vertical net magnetic 
flux strongly increases the level of random gas motion as 
well as the strength of the density waves in the midplane 
hence greatly reduces the maximum grain size achiev- 
able from grain growth, and may lead to collisional de- 
struction of planetesimals. In our laminar wind solu- 
tion, the disk midplane is essentially completely static. 
This provides the best environment for grain growth 
and settling, and is ideal for planetesimal formation and 
survival either via the streaming instabili ty mechanism 
(jJohansen et al.l 120091: iBai fc Stondl2010al) or the grav- 
itatio nal instability mechanism ( Lee et al1l2010l : lYoudinl 
l20Tll) . or via the vortices generated by the baroclinic in- 
stability (jLesur fc Papaloizo ul l20ll iLvra fc Klahdl20Tll) 
or the Rossby wa ve instability ( Lovelace et al.l Tl99£ : 
iMeheut et aTlf20Tl) . 

The efficient wind-driven accretion through the in- 
ner disk would change our understanding about the 
global evolution of PPDs. Most current models on the 
long-term evolution of PPDs adopt a phenomenologi- 
cal approach that treats the disk p hysics very roughly 
(|Zhu et al.ll2"0T0t iMartin et al.ll2012l ). In particular, un- 
der the framework of layered accretion in the inner disk, 
the presence of the dead zone often leads to inefficient ac- 
cretion and pile-up of mass. The pile-up gradually leads 
to gravitational instability which eventually drains most 
of the inner disk material onto the star, starting a new 
cycle. In the mean time, zero n e t ver tical flux global 
simulations by iDzvurkevich et all ()2010f ) found that the 
inner edge of the dead zone is likely to become a local 
pressure maxima which serves to trap solids. In real- 
ity, when more microphysics in the disk is taken into 
account, we see that as disk wind can be much more ef- 
flcient in driving accretion at the location of the conven- 
tional dead zones, mass pile-up and building up a local 
pressure maxima may be avoided. What determines the 
surface density profile of the disk is the spatial distri- 
bution of the magnetic flux, rather than the grain abun- 
dance, etc. Although our knowledge on the magnetic flux 
distribution in PPDs is very limited (jZhao et al.l 1201 ll 
iKrasnopolskv et aI]|2012D . and i t sensitively depends o n 
the internal structure of the disk ([Guilet fc Ogilviel2012|) , 
attention should be drawn to re-consider the structure 
and evolution of the PPDs focusing more on the trans- 
port of large-scale magnetic flelds using global simula- 
tions. 



6.3. Open Issues 

We have neglected the Hall effect in our calculation, 
which also plays an important role on the gas dynamics 
in PPDs. In the bottom left panel of Figure [SJ we also 
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show the profile of the Hall Elsasser number for our lam- 
inar wind solution. We see that the Hall effect at least 
dominates between z = 2H and 3H (within z = zL2H, 
a floor is added to Ohmic resistivity and AD, but not to 
the Hall term since it is not included in the computa- 
tion). At this location, wc have Ha ~ 0.1, /? '--^ 10, and 
the azimuthal magnetic field dominates the vertical field. 
We note that for the Hall MRI, the stability depends on 
the sign of vertical magnetic field, while for Ohmic re- 
sistivity and AD, stability is independent of the sign of 
Bz- We expect the magnetic configuration in our lami- 
nar wind solution to be stable to the Hall MRI when 
is negative, in which case we expect the magnetic con- 
figuration to be adjusted to account for the Hall effect. 
For positive Bz, our laminar solution may become un- 
stable to the Hall MRI a ccording to the calculations by 
iPandev fc Wardlj ()2012[ ). but the unstable wavelength 
may well exceed H due to the small value of plasma 
/3. Moreover, as vertical stratification is not included in 
their calculation, it is uncertain whether the disk would 
eventually relax to another laminar configuration or lead 
to sustained MRI turbulence, especially given the expe- 
rience in this work that an initially MRI unstable disk 
evolves into a laminar configuration in non-linear simu- 
lations. 

We have shown the the FUV ionization plays an 
important role in the wind launching process. In 
the mean time, the FUV p hotons are also important 
dirver of photoe vaporation (jGorti fc HollenbachI 120091 : 
lOwen et al.ll201^ . which is not considered in this work. 
Photoevaporation is considered as the main mechanism 
for dis k dispersal and largely determines the lifetime of 
PPDs ([Alexander et al.ir2006l) . We expect photoevapora- 
tion and magnetocentrifugal mechanisms to boost each 
other and enhance the mass loss rate compared with pure 
photoevaporation models, yet more realistic thermody- 
namics with proper treatment of heating and cooling 
need to be included, and eventually to make predictions 



to compare with observations ([Pascucci fc Sterzikl[200l 
ISacco et al.ll2012l ). 

The magneto-centrifugal wind is intrinsically a global 
phenomenon. As we have discussed before, uncertainties 
remains in our local shearing-box simulations regarding 
problems with the stability of the strong current layer, 
the mass loss rate, and the rate of wind-driven angu- 
lar momentum transport. Moreover, we have only stud- 
ied the gas dynamics for a MMSN disk at 1 AU, it is 
yet to conduct the same study at other disk locations 
to study how the wind properties depend on the radius. 
This will be presented in a companion paper, where we 
will show that the laminar wind scenario is likely to ex- 
tend to about 10 AU, while the MRI (in the presence of 
net vertical magnetic fiux) is likely to be the dominant 
mech anism for angular m omentum transport at the outer 
disk (H imon et aLll2013[ ). Meanwhile, another interest- 
ing problem arises on how the transition occurs from a 
pure laminar wind-driven accretion region to a turbulent 
MRI-driven accretion region. Global simulations are es- 
sential to resolve the problems and uncertainties men- 
tioned above, and are the ultimate way towards fully 
understanding the accretion processes in PPDs. 
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